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ABSTRACT 

A two-pool layer simulation model was constructed to 
describe solute transport through partially frozen, 
structured soil and a FORTRAN computer program written to 
implement the model. The model was based on a numerical 
solution of the continuity/chromatography differential 
equation. It can be applied to transport of non-reactive or 
adsorbed solutes through thawed or partially frozen soil. 

Effects of soil structure were described by defining 
two interactive solution pools, mobile and stagnant, based 
on the concepts presented by Addiscott (1977). Solute 
transport occurs only in the mobile pool and equilibration 
between pools is effected between infiltration events. 

Simulation experiments determined that the model was 
sensitive to stagnant:mobile pool size ratio, infiltration 
rates, adsorption coefficients and presence of a frozen 
layer, but not to layer thickness between 1.3 and ten cm. 
Simulations produced upwardly skewed solute profiles at 
higher solution flow rates and symmetrical profiles at lower 
flow rates. Frozen layers produced symmetrical profiles by 
reducing flow rates, thereby causing greater equilibration 
between pools. 

Simulations were compared with results of leaching 
experiments performed on a Gray Luvisolic soil during spring 
inftitration Of meltwater, susing “HOH “andy °*C-lindane as 
solutes. Effects of adsorption and frozen soil dominated 


over effects of soil structure on solute distribution. 
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I. INTRODUCTION AND LITERATURE REVIEW 


A. INTRODUCTION 

The transport of solutes through soil is a dynamic 
process involving physical, chemical, and biological 
phenomena. Estimation of degree and extent of leaching is 
important in studies in pollution, soil fertility, and soil 
genesis. 

One of the major challenges in the study of leaching 
through soil is determination of the quantitative 
relationships among the many factors which affect solute 
movement. In the field of chromatography the continuity 
equation has been applied with success to the quantitative 
description of solute transport through packed columns of 
homogeneous porous media (a list of symbols used is given in 
Appendix A): 

atdc/otvet+ Dbo(ea/dt) t= eulec7oz)etmDlo7c7#ozant 

This model of solute movement relates the rate of 
change in solute mass stored at depth, z, (a(dc/dt) + 
Db(da/at)) to solute flux density due to mass flow 
(-v(dc/dz)) and diffusion (D(d?c/dz?)). The rate of change 
in mass stored is divided into two components: mass 
dissolved in solution (c), and mass (a), held on the solid 
porous medium (by adsorption, electrostatic attraction, 
etc.). Thus the change in mass with time is expressed as the 


sum of partial derivatives. 
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Because soil is a porous medium that may adsorb solutes 
the soil column may be regarded as a chromatographic column 
in leaching studies. The model is useful because it can be 
used to assess the relative importance of the many factors 
which influence leaching. In the field, irregular patterns 
of leaching events and soil heterogeneity with respect to 
the parameters and variables used in the model make its 
analytical solution difficult. The chromatographic model 
can, however, be extended to provide a quantitative 
description of soil leaching, using numerical techniques and 
high speed computation. Further refinements can extend its 
use to predictive modelling. 

In short the chromatographic model is a powerful 
unifying concept in soil leaching studies because: i) it is 
a general description of solute movement; and, ii) it is 
quantitative. 

This study attempts to extend the model to solute 
transport studies in structured soils during spring thaw. 
Two major problems encountered under these conditions are 
the irregular "packing arrangement" of the soil medium and 
Spatial and temporal variation in patterns of water 
movement. The objectives of this project are: 

i) to construct a numerical model of solute transport 
through partially frozen structured soil, based on the 
theoretical chromatography equation; 

ii) write a computer program to implement the model; 


iii) test the sensitivity of the model to the parameters 


io 


easuloe diceby: yan. 2089 

nnufoo otdqsi~ezemoids & ze 6: 
ad ns> +i seusped Iujoru ee foton st at ne 
at0jos} yasm sit 3o sonsssoqubavizele # 
anisiiag zeipperys bleti pd? of vot 
of toeges1 dtie ytiensgotesed Lilet Sus as 
siam Isbom sfizv ni beawe asidatiev 8 


pei 
ishom >idaetootsmotds saT tint noisus ‘ee 


gHsb 


Wi 


stiinsup 6 sélivexg of Sebagixe saa 
7 


supindoess isoiremun gorau yentdoses ties Io'na 
7 
~43 Hneixe neo exaomentie®d s8hzTe .cSeitetug 


-prti{fLebom avitas 


i ishom sidgs:pesemorHs siz sr9ta'e 
‘ji ssauened eeibuse paidosef ffoe-nz tqeonee 4 
, bas. ,anem@svom asuloa to seitgizagse: 
» aa TS 

tehom att bretee of Biqnisitts yRuseeeiat 
oniiuh alioe Se rutzute ni asthuse Giege 
bnon seedt aebou Ss2sdans0one emeldoag no tent 
jfhem Lice sid te.."2nsmephstse prlaceg” 1siifpazyt 
to ahietied? n2 noltals ar isreqgus? Gis: 

‘ate spetoiwg sida )30, sevisostéo one 

sz 10 {sbom [astaeser & Souzten6s OF {at 
siz no bBeead fice Betusowsse assert yilsitsie¢ 
‘noltsEps yigeIpotsme mS Isoi 

iehom afd tnemelqmi of merpday Issuqmes = sak cua 


sr9t9msisq af? of Isbom sft 20 Yxtvisianes Sas see8 


and variables it uses; 

iv) compare the results of simulated experiments with 
results of field experiments to test whether the model 
Simulates solute transport in the field; and, 

v) determine what structural modifications are required 
and what further experiments are necessary to elucidate the 


important solute transport mechanisms which operate in soil. 


B. LITERATURE REVIEW 


1. Derivation of Chromatographic Equation 

The continuity equation was first applied to solute 
transport by chemists seeking to describe the behaviour of 
column separations (chromatography) (DeVault,1943; Kasten et 
al., 1952; Lapidus and Amundsen, 1952; Glueckauf, 1955). In 
1941 Martin and Synge introduced column partition 
chromatography (Zweig and Sherma, 1972), a concept which has 
proven useful in describing solute transport through soil. 
Chromatography can be defined as the study of movement of 
liquid mixtures through porous media. Soil 1S a porous 
medium and the soil solution is a liquid mixture; therefore, 
the soil system may be analogous to a chromatographic 
column. 

The classical theory of chromatography as it relates to 
Soils has been reviewed by Nielsen and Biggar (1962), 


Frissel and Poelstra (1967), and Kirkham and Powers (1972). 
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Leistra (1973) and Letey and Farmer (1974) reviewed 
pesticide transport models based on the chromatography 
equation. The former review stresses the mathematical 
techniques employed in computational models of solute 
transport, while the latter stresses the physics and 
chemistry described by the equation. Research related to 
Specific problems encountered in applying the classical 
theory to soils is the subject of this review. 

The following discussion develops the theory of 
chromatography as it applies to soil by deriving the 
continuity equation and reviewing contributions to the 
theory from both soil science literature and sources 


describing conditions paralleling those in soil. 


a. Mass Flow 

Mass flow, sometimes called convective flow, is the 
most important process influencing leaching of dissolved 
solutes and will be considered first in developing the 
chromatography equation. The continuity/chromatography 
equation, a statement of the law of conservation of mass 
applied to systems of liquid mixtures flowing through porous 
media, states mathematically that: any change in solute mass 
of solution passing through a volume element of porous 
medium must be accounted for by a change in mass stored 
within the volume element. This is illustrated in the figure 
I.1 and in the explanation following (adapted from Kirkham 


and Powers, 1972): 
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Figure I.1. Mass flow through a thin volume element (plate): 
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The mass, Ms, stored in the plate shown above, with 
thickness, Az; cross sectional area, A; and volume AAz, is 
given by: 

Ms = m(in) - m(out). 

Hence the change in mass stored with time is: 
AMs/At = A[m(in)-m(out)]/At 
= [m(z,tt+At)-m(z,t)-m(ztAz,ttAt)+m(z+Az,t) ]/At. 

The solute mass flux density, q(z), (rate of flow of 
mass through a unit cross sectional area of the plate) is 
given as: 

qdz)y sad dime \yidtdy/ A Caine. 

Therefore: 
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The mass flux out of the volume element can be 
expressed as: 

Qt zeAz) ee atzom tiidatz) 7daz Az. 
Therefore: 

OMs/ 0nm= Ae (zee G7) 4 GAZ <q (2) /0z')\) 

=A GANZ Ogee) /0'z 31)"5 

Sekutetmassetluxedensity,; gq, is related: to liquid 
Ssotutdond fluxedensity, v;sand Solute concentration, c: 

q=vec aie 
SO' 

OMs/ot = -(AAz)0(v-c)/dz S272 
Assuming that v 1S constant through Az, then: 

OMs/dt = -(AAz)v(ec/dz) Gos. 

Ms can also be related to solute concentration: 

Ms = V°c, 


where V = Vp (i.e.Volume = pore volume) 


a*Vb (porosity X bulk volume) 
=NaNUA AZ. 
Som tha ts 

dMS/dt = a(AAz)dc/dt, 
if a does not change with time. 

Therefore if solute is stored in solution only 
(substituting the above expression for the left side of 
equation (6)): 

(RAZ ooc/ Ob) = 1 - NAAZ) UV. COC 70S) hn Or 

aloc/ot) =] > uNoe/ oz) G2 


which equation describes solute transport by mass flow only 
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and solute storage in solution only. The other main 
assumption is that the porous medium is homogeneous with 
respect to flux density, porosity and moisture content. Note 
that v in this case is the liquid flux density, or Darcian 
velocity, and not the pore velocity, v, as used in the 
equation given by Kirkham and Powers(1972). To convert to 
their equation we need only note the relationship between 
pore velocity, v, and Darcian velocity, v: 

y 2Ovpa)¥So% 
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b. Diffusion 

SOLUute fluxsaliso occurs as a result of diffusion in 
response to concentration gradients. When a band of solute 
is applied to the soil surface, as when a pesticide is 
sprayed onto soil, a concentration gradient 1S set up and 
diffusion through the soil solution is induced. Fick's Law 
States that the diffusive flux, qd, is given by: 

Cde=t- DCG /dzs 
where D is the diffusion coefficient of proportionality. 

Diffusive and convective fluxes can be summed to 


express total flux,’ gt: 


Gt =40C 770d 


qce7- -DpVde7azo. 

By substituting qt for q equation (3) becomes: 
8Ms/dt=n-(AAz)dlqe - Dldc/daz) ]/az 

=o-(AAz)idge/azit CAAz)DGd2¢/0z2), 
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Recalling that convective flux, qc, is given by 
equation (4) (qe = v-c) and provided v is constant, then: 

a (OG / ati) f=m= “v Gac4ezin+iD (6307827) Gag pe 

Kirkham and Powers (1972) and Leistra (1973) present 
the form: 

d0/o0ti= i-BrvGic 7dz) enebt (02¢/bz34), 
whichsis ’identicaleto"(8)cifier =lv/arvandvD'= °D/a: 

The factors which influence D are given by Leistra 
(1973) to be volumetric moisture content (8), geometry of 
the porous medium (tortuosity, etc), and partition of solute 
between solid and liquid phases, as well as the value of D 
in water (D.). This study treats partition separately so D 
will be not be regarded as a function of partition 
Deapanecer se 

Frissel epeal’. 101970) define Deas the ®product of °D,, 8, 
and tortuosity, but the theoretical basis of this expression 
was not given. Shearer et a]. (1973) define D as: 

BD i=ni6/a)40* "Dp, 

The latter expression does not require determination of a 
tontuosutytfiacton. 

The preceeding discussion applies only to longitudinal 
diffusion. Transverse, or intra-aggregate (intrapedal) 
diffusion will be treated later under the effects of soil 
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c. Dispersion 
Spreading of the solute band is also caused by 


"hydrodynamic dispersion". While the term diffusion refers 
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to the spreading of solute particles within the soil 
solution, dispersion refers to the longitudinal spreading 
caused by spatial variation in fluid flow rates. 

Although the chromatographic equation deals with an 
average flow rate, solutions do not move at uniform rates 
through soil. Even in completely homogeneous media some 
mixing of fluid occurs because of adhesive and cohesive 


forces along pore walls(Figure 1.2). 


Figure I.2. Hydrodynamic dispersion due to "drag" ina 


capillary tube. 
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(Arrows represent the magnitude of fluid velocities and 
distance travelled by fluid "layers".) 


This diagram represents "laminar flow" resulting from 
shearing of fluid in the longitudinal direction. This occurs 
because the attraction of the solid for the liquid retards 
movement close to the walls. Cohesion between liquid 
molecules resists flow away from the pore walls. Maximum 
flow rates occur at the greatest distance from the pore 
wall. 

In the capillary bundle model (Nielsen and Bigger, 


1962: Scheidegger, 1974) this is the mechanism of dispersion 


“tioaveds 7 
—— Lenibutt 
.zede2 wolt 

ns fitiw sissh- nobisups 3 bt 3.50" 
aeis1 migiinu ts avom Jon ob tativion ¢ 
smoze 6iben avosnepomed yisselqmoo n 
syizgeden haa avieed&s io saussed a1wsame 

a oe situpitieliaw 

aan 


7 
s ni "gs1b". oF -eub nolaveqeib. pimanybowbglt « 


: ie A Pee 

gsstzisolsv Giuli to abu singes saz 3a | wo’ 
(,"ayeyei" Bialt ed Pa yn ie ete 

mo12 enisiveas2 ‘wold wenimel* 2insGs3Gs7 astgeih -eidt : 
syo90 etd? .noticetib. laaPbuzipnol eda ni Sines prt 
hiupil eds s02 BStioe ads 10 norte7sae oA 
Siupil naewted netasdo> .ellew ed’ 03 SmQEoe 
numixeM .ellew exoq edd most ysus woll etefaege 


etoqg sd3i moti sones2ib testae10 ody Js W590 


a! 


teppid bas nseleit) Lebom sfibaud val Liges edz 
: vs ra ; 


norateqeid ia aa Enedoad edz ai aids (ste! 2 >: 
7 


ra 


10 


if flow is laminar. The degree of dispersive mixing in such 
a model depends on the pore radius and fluid velocity. In 
SOil, however, the pore geometry is much more complex so 
pore radii cannot be well characterized by an average value. 
Therefore, characterizing liquid fluxes by an average value 
is frequently not realistic. 

Amouzegar-Fard et a]. (1982) tested the effects of 
probabilistic versus deterministic approaches to fluid 
velocity. They used normally distributed values of moisture 
flux to calculate 2000 separate solute concentration 
profiles in simulation experiments. They found that the 
average of the 2000 profiles differed sharply from the 
profile calculated using the mean value of moisture flux. 
This experiment provides a good illustration of the effects 
of dispersion on a larger scale than the capillary tube 
model. 

Dispersion is treated in the classical equation ina 
manner analogous to diffusion (Kirkham and Powers, 1972): 

G@arsp) Maas (ac7072)), 
where Bb iS the coefficient of proportionality “(cm-/d). 
Because both diffusion and dispersion result in a spreading 
of solute they are commonly combined so that: 

Giausp ye q(t rtemea—e-Dsprtdce/dz) « 

Lumping D and E together as either Dspr or Da (Da = 
apparent diffusion or dispersion coefficient) necessitates a 
minor change in equation (8): 


aloc/dt) = -vldc/dz)+ Dspr.(o7c/o0z7) (9). 
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Diffusion perpendicular to the direction of flow 
reduces the effects of dispersion. Consequently dispersive 
and diffusive effects are not additive, making resolution of 
Dspr into D and E difficult. Nevertheless, Frissel et al. 
(1970) experimentally determined D and E, and defined Dspr 
as: 

Dspr = v-E + @°T-D Oo; 


where T = tortuosity (dimensionless), 


8 volumetric moisture content (ml/cm?). 
OQtlalitatively,°(10) states’ that diffusion effects are 
more important for low flux densities, while disperion 


effects are paramount for high fluxes. 


d. Partition: Adsorption-Desorption 

Up to this point phenomena which result in movement of 
solutes have been discussed. One of the main factors which 
cause retention or "holdback" of solute 1s adsorption. A 
Substance is Said to be adsorbed if the concentration of the 
substance (sorbate) in a boundary region (solid surface) is 
higher than the interior of the contiguous phase (soil 
Solution) (Tinsley, 1979). 

Chromatography utilizes the adsorptive partition of 
solute between mobile and stationary phases to separate 
Solutes within or from solutions (Johnson, 1972). 
Separations occur in a similar fashion in soils. Soil 
Scientists were the first to use isotherms to represent the 


adsorption-desorption partition (Giles, 1970). 
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For adsorption of many pesticides by soil the 
Freundlich isotherm has proven sufficiently accurate (Kay 
enaemilickw 1967 mmodgdsonmenr al... 1970- Lindstrom and 
Boersma, 1970). It describes an empirical relationship 
between amounts of adsorbed and dissolved substance at 
equilibrium, in which: 


a=" Ka0ds+cac) is, 


= 
Sr 
0) 
ual 
49) 
0) 
i 


amt of adsorbed substance (ug/g soil), 
c¢ = solution concentration (ug/ml), 
ke=—propontronality constant (ml soln/gq soil), 
1/n 


The exponent in this equation has been found to lie between 


empirical exponent. 


0.70 and 1.00 for pesticide adsorption in most soil systems 
(oeretray 11973). 

Giles et a]. (1960) classified Picaconn isotherms on 
the basis of their shape, dividing isotherm shapes into four 
classes and five subclasses. Freundlich iostherms fall into 
"class L" (Figure I.3). They hypothesized that for this type 
of isotherm solute particles are most likely adsorbed flat 


against the surface of the sorbent as in figure 1.4. 
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Figure 1.3. Freundlich (Giles Class L) isotherm. 


a =akadse ce. 


equilibrium soln concen (c) 


Figure 1.4. Schematic diagram of adsorption for "class L" 
(Freundlich) isotherm. 


adsorbed substance (sorbate) 


a ae oe ee ee 


particle surface (sorbent) 


At lower concentrations the plot of a versus c in 
figure 1.3 is almost linear. Consequently for sparingly 
soluble pesticides a linear isotherm is commonly used. 

To accomodate adsorption-desorption in the classical 
equation recall (6): 

OMs/dt = -(AAz)v(dc/dz). 

If storage of solute occurs in the adsorbed phase as well as 
in solution then: 

Ms = a(AAz)c + Db(AAz)a, 
where Db = dry bulk density of soil (g/cm*) and all other 


symbols as previously defined. 
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Substituting this expression into (6): 
aldc/dt) + Db(da/dt)= -v(dc/az) (Gipals 
This can be converted to the equation of DeVault(1943): 
aickde/oV:) = 40GN/ON = -0c /oz 
wier 6 jun= di /A GOV / 0b). eormiot = 0aVi/ uA) > 
Ovec= aDbi> Ablavesvand 
a' = a*A (porosity per unit column length). 
DeVault has used dV aS a convenient measure of time. This is 
valid when v is constant. DeVault (1943) used the general 
isotherm: 
OF i=yni ee ersothat’: 
lens’ iecd £Cci)sfociioc/oV.= -dc faz. 
Spreading by diffusion can also be included in (12): 
a(ac/at) + Db(da/at) = -v(dc/dz) + D(d2c/az?) Ga 
Dividing through by a and allowing for the different 
dimensions of a used, one obtains an equation identical to 
that of Lapidus and Amundsen (1952). 


Relating a to c by a linear isotherm (a = Kads-c) (13) 


becomes: 
akec/oL)s + Dboekaasuce/ dt — ude /dz) steDUasc (oz), OY, 
Cao+eDb«Kadsiocyon = =i loc/oz) + Dld2c/o0z-) (a8): 


Again dividing by a to change dimensions (14) can be altered 
to the equation used by Kay and Elrick (1967). Equation (14) 
has been developed for saturated flow, that is, 6 = 6S = a. 
For unsaturated flow 6 should be substituted for a. 
Substituting Dspr for D, as in (9), the chromatographic 


equation becomes: 
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(6 + Db:Kads)dc/dt = -v(dc/dz) + Dspr(d?c/dz?) C15 )*. 

So far the kinetics of the adsorption-desorption 
reaction have been ignored. It has been assumed that rates 
of adsorption and desorption are high enough, relative to 
the liquid flux, to permit transverse equilibrium between 
adsorbed and dissolved solute. This may not always be a 
valid assumption. Lapidus and Amundsen (1952) dealt with 
rates of mass transfer between phases, and Kasten et al]. 
(1952) described the case in which rates of diffusion to 
adsorption sites within solid particles limited 
adsorption-desorption, which may be a more realistic model 
to describe leaching in soils. Rao et a]. (1980), Addiscott 
(1981) and Nkedi-Kizza et a]. (1982) have developed 
numerical models describing the effects of finite 


intra-aggregate diffusion rates on solute transport. 


e. Sinks and Sources 

The continuity equation can be written to include sink 
Or source terms: 

Q(a + Oc) /at = -a(vc)/az +°S (16), 
where S = sink or source term due to salt precipitation or 
dissolution (Magdoff and Bresler, 1973). 

Rolland and Frissel (1974) used a first-order decay 
reaction rate law to describe pesticide decomposition in 
soils: 

depot i=an- kc Cat), 
where k is the decay rate constant (1/d). 


This is an approximation of the Michaelis-Menton equation 
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for very low concentration, as is uSually the case for 
pesticides in soil. They have also used a first order decay 
rate law to describe irreversible adsorption of pesticides. 

Equation (17) can be solved for c(t+At): 

ChttAt) {ac Ge jexol-heAt:) 
Therefore a simple way to deal with decay of pesticides is 
to multiply the solution obtained from the chromatographic 
model by the factor exp(-k-At) at the end of each time 
period,. At. 

The non-equilibrium case for adsorption-desportion was 


described in a Similar manner by Lapidus and Amundsen 


(1952): 

da/az = K,c — ‘kea (18). 
where kK, = rate constant for adsorption, 

kz = rate constant for desorption. 


£. Soil Heterogeneity 

Effects of soil heterogeneity on leaching were 
mentioned in the subsections on diffusion and dispersion. 
The most obvious feature of heterogeneity is the vertical 
Stratification which characterizes all soils. This 
Stratification results in variation in the physical and 
chemical parameters used in the chromatographic model. Many 
of the parameters were assumed to be constant during 
derivation of the model described in previous subsections. 
The situation in soil is better described by the following 


general model (Leistra,1973): 
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oMS/0t = —-o(v-c)/oz+ 0(D+dc/0z) /dz 

+ 0(E-dc/dz)/az - (8M/at)cons (19), 
where (08M/at)cons = consumptive term which describes affects 
of assimilation, decomposition, and fixation; 

OMs/dt = 0(6°c)/dt +d3(Db-a)/at; 

ange where 6, Db, iu, D,. and E can vary with depth. 
Solving this equation is one of the current challenges to 
research in soil physics. 

Vertical Stratification is not the only aspect of soil 
heterogeneity of concern, however. Structural heterogeneity 
(pore distribution) influences diffusion and dispersion. 
Channelling of solution through macropores which represent 
only a fraction of the total porosity, results in greater 
depth of solute penetration than is estimated on the basis 
of total porosity (Tyler and Thomas, 1981). 

The soil solution may also exhibit variability. Smiles 
and Gardiner (1982) have dealt with heterogeneity of the 
soil solution as it affects anion transport. They divided 
Ehewsoid SOluULLON Into stwO.pools*+ one: accessible sand, one 
inaccessible to anions. In this way they describe anion 
exclusion in the solution very near to the surface of 
negatively charged soil particles; inaccessible water was 
defined as a layer 9A° thick. Their model appears as: 

O(Ga-ca)/ot = d(B<-dca/dz)/oz — d(v-ca)/oz CS0)2 
where the subscript "a" refers to values for accessible 


water only. 
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Addiscott (1977) has used a similar approach to treat 
anion exclusion in structured soil. The case where a solute 
is excluded from a region is a simple illustration of 
heterogeneity of the soil solution. Although this study does 
not consider anion transport, the idea that the soil 
solution can be divided into separate phases, or pools, in 
order to treat the effects of soil strtcture on solute 
transport is central to the model which is developed herein. 

Where intra-aggregate diffusion rates exert a kinetic 
influence on adsorption-desorption (as in the case of 
relatively high v) radial diffusion within aggregates must 
be considered separately from inter-aggregate rates. Kasten 
et al. (1952) considered this phenomenon. They solved an 
equation parallel to the classical model of vertical solute 
transport for radial (intra-particle) transport. Addiscott 
(1977) has simplified the treatment of diffusion rates into 
and out of aggregates in his two-phase model of the soil 
solution. He assumes that during infiltration events there 
is no (diffusional) mass transfer between mobile 
(inter-aggregate) and stagnant (intra-aggegate) solution 
phases. Between infiltration events, however, complete 
transverse equilibration of solute between phases occurs. 
These simplifying assumptions avoid the problems of 
determining aggregate sizes and diffusivities (D) which vary 
with Gepth and; «in the case of D, with moisture content. It 
is necessary, however, to examine the validity of these 


assumptions, since in the range of intermediate flow both 
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may be invalid. 

Skopp et a]. (1981) have used an approach similar to 
Addiscott (1977), in that they defined two regions of soil 
solution but regard both as mobile. Each solution pool moves ~ 
at a different rate and a small degree of solute mass 
transfer occurs between them. Their model would best 
described cases where soil solution flows around and through 
Soil aggregates, while that of Addiscott (1977) describes 
cases where solution flows around aggregates. 

The chromatography equation can be used to describe the 
effects of mass flow, dispersion, diffusion, adsorption, and 
decomposition on solute transport. For solute transport in 
Soil, heterogeneity resulting from stratification, and pore 
geometry (soil structure) require modifications of the basic 
equation. The information and assumptions required for 
solution of various forms of the model are discussed in the 


next section. 


2. Solutions to the Chromatographic Model 


a. Analytical solutions 

The techniques of differential and integral calculus 
can be employed to solve some forms of the chromatographic 
equation. These solutions are exact and usually simple to 
use. Use of analytical methods to solve the model, however, 
is severely limited by the complexity of leaching problems 


encountered in the field. 
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For instance equation (19), which is general, is too 
complex for analytical techniques. Most analytical solutions 
assume that v, 6, and Dspr are constant so that (19) is 
greatly simplified. Further simplifications concerning 
boundary conditions and feed functions are necessary for 
application Cffanavytacal) solutions (Leistrayi973) In 
light of the complex situations which occur in soil these 
Simplifying assumptions cannot be justified and numerical 
Solutions must be sought. 

Analytical solutions are nevertheless instructive. They 
can be used to analyze the effects of variation of 
parameters (e.g., Dspr, partition coefficient), soil 
physical properties (6, Db, a), and feed variables (c(t), 
v). They are also useful in assessing the error introduced 
by numerical methods. This can be done by using the same 
values for all variables and parameters in solutions which 
are available in both numerical and analytical forms. 
Concentration profiles can be calculated from analytical 
solutions with the assistance of only a pocket calculator 
and mathematical tables, whereas numerical solutions require 
computer assistance to accomplish the volume of computation 
involved. 

Kirkham and Powers (1972) derive solutions for some 
simple cases described by the classical model, such as the 
general solution for a one-dimensional dispersion model 
(dc/dtac=Gbl0%e/d277) )e This equation is solved fon 7a moving 


coordinate system (t,,z,) and then transformed to the fixed 
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coordinate system (t=t,; z=z,+vt), that is, the dispersion 
model is superimposed on the piston-flow model of mass flow. 
They apply this solution to two elementary boundary 
conditions:(i) dispersion of a displacing front; (11:1) 
dispersion of a "Slug" (narrow band) of solution. They also 
provide the Carslaw-Jaeger solution for the mass flow - 
auitusvond equation (0e7 0thizerilec/ozjet DldFcC/0z 24 )r 

Both D and E can be calculated from experimental 
breakthrough curves using analytical solutions (Kirkham and 
Powers, 1972). The solutions mentioned above illustrate the 
differential and integral calculus required to solve the 
model analytically. Neither of these equations considers 
partition between solution and adsorbed phases and therefore 
applies only for non-reactive solutes, such as chloride ion 
or tritiated water. 

DeVault (1943) provides an analytical solution for the 
model describing mass flow with adsorption-desorption 
(non-linear isotherm). Instantaneous transverse 
adsorption-desorption equilibrium is assumed and spreading 
by diffusion or dispersion is neglected. The solution would 
be valid only for a range of fluid velocities for which 
diffusive spreading is insignificant in the lower limit of p 
and dispersive spreading is insignificant in the upper limit 
of vp. Fluid velocities must also be low enough to exert no 
kinetic influence over adsorption-desorption. 

Lapidus and Amundsen (1952) consider an apparent 


contradiction in DeVault's assumptions, namely, no diffusion 
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and instantaneous adsorption-desorption equilibrium. They 
reasoned that when mass flow rates were low enough for 
point-wise equilibrium at all levels in the column, then 
longitudinal diffusion rates would be significant. 
Similarly, at higher flow rates which eliminate the effects 
of diffusion, point-wise equilibrium in the column is a 
questionable assumption. Therefore they solved the model: 

dc Gt /PEl ll falvdnfott= =v (oc/dz) + Dtarc/az 4) 
where dn/dt is the rate of exchange of mass between solution 
and adsorbed phases. They solved this equation for 
equilibrium and non-equilibrium conditions. 

Kasten et a]. (1952) considered the case in which 
finite intraparticle diffusion rates limited 
adsorption-desorption. They defined an/dt in terms of the 
rate of intraparticle radial diffusion. 

All of the above analytical solutions apply to boundary 
conditions and solution feed functions which can be 
carefully controlled in vitro. Research in soil under field 
conditions allows no such luxury. Moisture flux under field 
conditions is a variable function of time while severe 
discontinuities may occur in soil properties. The analytical 
techniques used in the above solutions are based on 
continuous functions. Any rigorous treatment of solute 
transport through soils must be capable of treating the 


complex variation which occurs in the field. 
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b. Numerical Solutions 

Numerical solutions permit consideration of more 
complex soil properties and leaching conditions than can be 
treated by analytical techniques. The first step in most 
numerical approaches is to divide the soil column into 
discrete layers so that a difference equation can be 
written: 

DbC0a/dtpe +0 ocy0t) e=ag-vAcyAze+ DsprAec¥Aze (2h). 
For non-reactive solutes: 

Olde /dtym-favAcyAznsrDsprA2cyAz2. 
If Azgeiseset = 1, then: 

Gep=a(170) (jvAcetiDsprA2zc) at: ox¢ 

tene Poctaad ees A046 hav dGeteDSprAécids 
whéererAce=eects zien etsy zit); 
and Ate = /oute zed) ga /2echkt+2) Or /e(t , z-de 
So thaté 
Cibhinz) kaec Gt ewe: FOn (4 / beter oe) ec olty 2) ldt PAM 

Frissel et a]. (1970) have used this approach, solving 
the integral using a 4th order Runge-Kutta formula. This 
formula is very accurate but is based on the presence of 
very smooth concentration versus time curves (feed 
functions). For discontinuous feed functions other numerical 
integration techniques can be used to solve (22), of which 
rectangular integration, though less accurate, is the 
Simplest: 
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Leistra (1973) developed a finite difference equation 
based on rectangular integration for the classical model 
describing leaching of non-reactive solutes: 

Clt+1, 2)0= eo hazyterAt/Agio-cie azet)--peicty ez) 

SPOSDERCGE TZ Ie se (t 2-1) )/Az 
TDS OUT CW ce meee C(t rz) ize 2au 
In agscGlutiomoorzthismtype initial conditions define c(t,,z) 
Portal bizeuNnextecet oraz) hestideterminedatromvyc (tgipz) tat ,edch 
layer and v from equation (23). This process is repeated for 
each time step for the period of time for which the solution 
is desired. 

Variations in soil properties and flow rate are often 
important in leaching studies. Recall the essential 
differences between equations (15) and (19): 

GG foeDbakads )oc/dt =l=au(deZoz)il+ Dspr( 440/324) (15) 

0(O@°c)/dt + d(Db*Kads-c)/dt = 

-d(v-c)/az + O(D°dc/dz) /dz 
+ 0(E°dc/daz)/az + (38M/dt)cons €4:9:).. 
Equation (15) is specific for columns in which 86, Db, Kads, 
v, and Dspr are constants and can be handled with analytical 
methods. The general equation (19) describes vertical 
heterogeneity with regard to these properties, but requires 
numerical methods to yield a solution. 

Leistra (1973) solves the general equation for leaching 
of non-reactive solutes: 

a6 <x) Zatle= 8- Or ce) 7idzmet to (DsprsdeA0z)i/ dain koarnse0-c} 


where kcons:@:c is a decay sink term. Four coefficients (6, 
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vy, Dspr, and kcons) can vary with depth. 


The finite difference solution is: 


Cte ez) eae G aie 21) /ORZ ic (te 21) 
Me tee Alec ernie hae axle) AOC zi) users (a / Oz) VedteZ)) 
zy Cue Ct potsy:, 
Where i, =CvAt/ Az. 
RAUz) = sDSpriézy Atyaz:. 


c(t+1,z) can therefore be determined at any time, t+1, from 
ett athe concentration profile at time, tt, and ‘the 
necessary soil properties defined for each horizontal layer, 
Zeon ncne SOlDMeCoLumn 

The numerical solution of Frissel et aj. (1970) treats 
variable 6 and Dspr, but deals only with constant water 
flux. Addiscott's (1977) model of a structured soil system 
allows for variation in flux but only very ned variation 
in 6. Burns' (1974) model solves the solute transport 
problem for variable ®@ and flux. Dispersion and diffusion in 
the models of Burns and Addiscott are effected by numerical 
spreading only. Leistra's (1973) model deals with vertical 
Stratiticacion out not swith ysomlestructure. 

Published information about simulation models is 
frequently insufficient to assess the numerical techniques 
or to reproduce the simulation experiments. De Wit and van 
Keulen (1972) and Frissel and Reiniger (1974) published 
detailed manuals and programs describing a variety of 
Simulation models of transport processes in soil written in 


CSMP. The former include solutions of the moisture flow 
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distribution, while the latter require that the flux be 


known sfromecother ‘solutions or can be controlled. 


c. Plate Models and Rate Models 

Most of the models mentioned above are based on 
rigorous physical laws. Some assume that there are no rate 
limiting steps (DeVault, 1943). Others include specific rate 
laws (Kasten et a]. 1952; Lapidus and Amundsen, 1952) to 
describe the effects of diffusion, dispersion, finite 
adsorption-desorption rates and chemical reaction rates. 
Those models which contain rate laws to describe kinetic 
effects are called "rate theories" by Frissel and Poelstra 
(1967). 

There is another important chromatographic model which 
relies on the theoretical number of plates ina 
chromatographic column, or the "height equivalent to a 
theoretical plate" (HETP), to determine the solute 
distribution. Frissel and Poelstra (1967) call models of 
this type "plate theories". The number of theoretical plates 
is equivalent to the number of theoretical equilibrations, 
"n",. This value can be used to calculate the solute 
distribution in the column as follows (dJohnson’, 1972): 

C= ninl/ Lina ab B SY los ee (B/C Beale (24), 


= concentratzon in the 1th plate, 
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effective partition coefficient (Db-k/6), 


to 
Wl 


number of theoretical plates or equilibrations. 


= 
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The above expression is the ith term in the binomial 


expansion of [1/(B+1) + B/(B+1)]", where n can be calculated 
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as follows (Johnson, 1972): 


m= Sp(L“p) /w7, 


where p = position of the solute peak concentration, 


L = length of column (depth of penetration of leaching 
solution), 
w = width of solute band at point where the solute 


concentration is 36.8% (1/e) of the maximum concentration. 

Tmmtthiiss conte stiing cor HETPseni/ntias asedia's 2 
parameter of the leaching system. Glueckauf (1955) developed 
the plate theory for continuous feed functions (liquid flux) 
and showed that HETP is a measure of the spreading 
coefficient. A theoretical value of Az which when used in 
the models of Johnson (1972) and Glueckauf (1955) accurately 
predicts spreading of the solute band ee particular set 
of conditions. 

For finite difference models, such as those of Leistra 
(1973), the value chosen for Az also influences the 
calculated spreading of the solute band. In finite 
difference models, where concentrations are uniform within 
each layer, spreading will be greater when calculated on the 
basis of larger Az. This effect is known as "pseudo 
dispersion" (Leistra, 1973) or "numerical dispersion" 
(Frissel and Reiniger, 1974). Numerical dispersion can be 
reduced by using relatively small values of Az, which 
technique requires more computations. Alternatively Az can 
be chosen so that the numerical dispersion produced by the 


Simulation model corresponds to actual spreading produced in 
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leaching experiments. This value of Az can is determined 
from the results of a leaching experiment. Using Az values 
in place of rate laws to describe spreading is invocation of 
the plate theory. It is sometimes difficult to apply plate 
models to soil because the HETP requires approximately 
constant values of flow rates and spreading coefficients. 

It 1s also possible to prevent numerical spreading by 
selecting the time step of the Simulation so that the solute 
concentration profile is shifted down exactly one layer 
thickness, Az, at each time step (Leistra, 1973). This 
allows the modeller to introduce an accurately known amount 
of spreading. For field situations this method would require 
a variable time increment. 

Values of Az used in simulation models may also exert a 
pseudo-kinetic effect on simulation experiments. The case 
where flow rates exert a kinetic effect on partition of 
solute has already been discussed (Lapidus and Amundsen, 
1952). In the layer model of Addiscott (1977) this effect is 
Simulated by using two different leaching routines. A slow 
leaching routine (SLR) allows local equilibrium within each 
layer of thickness Az. A fast leaching routine (FLR), which 
Simulates intense rainfall events, limits the interaction 
between mobile and stagnant solution pools, thus simulating 
the kinetic effects of flow rates on partition between 
solute pools. Entry into the fast leaching routine depends 
on the quantitative relationship between Az and the rainfall 


input, AV (which depends on At). Therefore At must be 
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compatible with Az for the simulation model to correctly 
predict depth of leaching. 

Normally Az is determined by a leaching experiment. One 
requirement that Addiscott (1977) sets for simulation models 
of leaching, however, is that no prior leaching experiment 
be required. He gives no @ priori criteria for his choice of 
Az. He has, however, determined the Az for the best fit 
between calculated and experimental values, thus, in effect 
performing a leaching experiment to test his choice of Az. 

In a later version of the model Addiscott (1981) 
included a specific rate law for diffusion between mobile 
and stagnant solute pools. Thus he transformed the model 
into a rate model, and precluded the need for a prior 
leaching experiment to determine Az. The specific rate law 
can define mobile-stagnant solute partition for any Az and 
At (or AV). 

The above discussion illustrates the effects Az can 
have on simulation models. To avoid numerical dispersion, Az 
must be relatively small. HETP can be used in place of a 
rate law where leaching can be characterized by an average 
HETP. When leaching is slow relative to partition rates, Az 
need be chosen with regard to its effect on numerical 
dispersion only. 

Addiscott (1977) assessed the sensitivity of his model 
ofrleachingoathroughtstructuredisorlsutocAz waluesniroma3.25 
to 13 cm using time steps of one day and concluded that 


little benefit was gained by reducing Az below 3.25 cm. 


e¢ 7 


= pie ay 


yljoe1102 03 Lebom 5 


Nai ; 
ae uss te a 
siabom noicgsicmie a2 @292- (t | st08 ibe 


tnemitzeqxs en nidesekl w38lTy cme ‘38 vod 


onO .tremizeqxe peidoss! so yeoeeme 


estod> eid “st stiediz inate S on zavipos ah 


aN 


tii»s2sa e139 101 zi ae) banterrsieb ] ra 


inatia ni ud? ,esdlav iascatmies one, 


spies aid se93 ot s0amizequs enidoaal . 8 


sAio ri 
('82l) ttovgibSA febom aie i@ sores 1 ee 


se 10? Bean adi bsburlosta@-ane’, 


inf 
mt 
| 
- 


wal ste1 tibizseqs aa? 2 ociotese inemixeque 


bas s& vns tot notaitieg stutoa tnongate-eiidow ai 3 


a 
s 
wu 
, 
wW 
a 
et) 


serseuili sotezuseib ayeda eee 


fortetecgaib issizanya Bices of -eiedom not telumes 


- f * Ss _ * o aad = A & 
s i0 ansiq nt beau ad nso STSH .tleme J svigeesg So am 


A 
speaTavs os yd besisasgoe tsdo 30 585 onisosel siadiw wells 
3 7 
“4 .ee¢e% noisittse of svigsiey wole ef patioam® 


ieaitemua so ¢ae2320 ath of Srepées délw ageetey 
. (ine ae? 


cabom gid 20 vtivitienss sag Beazeces (YVET? t4o0a 1 bBA 


> £ mot? esulsv sti o@ aliog beiwsousse fioupiAd prid 
: : x 
isdi bsdulones ans sb ‘soe. to eqeda emit Dake me 


mz 28.2 woled sh pricuiem ya bentep. 2ey J 


; 
an 


30 


Thoma et a]. (1978) determined that their plate model for 
tritium leaching through sand dunes was insensitive to HETP 
variation for HETP values below ten cm. Burns (1974) assumed 
instantaneous equilibration between stagnant and mobile 
solute pools. A Az value of 3.75 cm. yielded simulation 
results which were in good agreement with experimental 
patterns of chloride leaching. 

Thornthwaite et a/J. (1960) used the analytical plate 
model presented above (24) to determine radiostrontium 
leaching patterns in laboratory columns of four texturally 
different soils. They used an HETP value of 1.27 cm. for all 
four soils. Their predicted profiles showed a good 
correlation with measured concentration profiles. 

Numerical solutions of the chromatographic equation are 
better suited to solute transport in soil than analytical 
solutions because they can be applied to the complex 
variation which occurs in the field. Use of numerical 
solutions requires special attention to the error 
introduced. Error can arise from both the time step and 
layer thickness and sensitivity of numerical models must be 


assessed for both sources. 


3. The Transport Model Applied to Structured Soils 

The structural organization (aggregation) of soil 
particles controls pore size distribution which in turn 
influences flux distribution. The Hagen-Poiseuille equation 


States that volume flux through a pore is proportional to 
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the fourth power of the pore radius (Scheidegger, 1974). 
This relationship between flux and pore radii indicates that 
water contained in small pores may remain nearly stationary 
while water nearby 1S moving rapidly. Some of the water 
present in the field, therefore, may not participate 
Girecehy intleaching (Addiscott)e1977)< 

The pore distribution of strongly aggregated soils may 
be defined by two distinct size ranges: intra-aggregate 
pores, or micropores, and inter-aggregate pores, or 
macropores. Zimmerman et a]. (1965) approached soil 
structure aS a Simplified dispersion problem; that is, they 
characterized the flux distribution by two fluxes, v,=0, and 
v2>0. The weighted average of v,;, and vz would be the mean 
volume flux. 

Addiscott (1977, 1978) treated solute transport in 
structured soils by defining two solution phases on the 
basis of their mobility: mobile, or inter-aggregate 
solution; and retained, or intra-aggregate solution. A 
fraction of the solute in each layer is considered to move 
at each infiltration event and the remainder to be retained 
within its respective layer. After each downward 
displacement of solution an equilibrium is established 
within each layer by equalizing the concentration between 
pools. The assumption of complete equilibration between 
intra- and inter-aggregate solution may not be justified in 
all cases, however. In a later model Addiscott (1981) 


defines interaction between moble and retained pools by 
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incorporating a diffusion rate law (Fick's law) into the 
model. Solute exchange between pools is controlled by 
diffusivity, aggregate size, and intra-aggregate 
concentration gradients. Rao et a]. (1980) dealt with soil 
micropore regions as a diffusion sink/source so that 
movement into or out of aggregates is defined by a 
Sink/source term in the classical chromatographic equation, 
which they solve numerically. The sink/source term is 
defined by radial diffusion rates. 

The problem in using the two pool concept is 
development of a criterion to define the pore size range for 
each pool, or the volumes assigned to each pool. Addiscott 
(1977) divides the pools at the 200 kPa soil moisture 
tension, with the upper limit of the mobile water pool set 
at fed capacity (5 kPa) and the lower limit of retained 
water at the volumetric midpoint between wilting point (1500 
kPa) and dryness. Thus he implies the existence of a third 
pool which does not interact with the other two. His 
reasoning is that this fraction of Boia water 1S associated 
with extremely fine pores and diffusion within these pores 
is slow. Since he was modelling anion transport he assumed 
that this pool excluded the solute. 

Addiscott (1977) presents the 200 kPa division as 
st muctlyyempir icalijaithoughs he istatesythatoathishouldt bes the 
same for all soils, but may vary with the solute. He does 
not indicate why he has used five kPa moisture tension for 


field capacity, rather than the conventional ten or 33 kPa. 
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Nkedi-Kizza et a]. (1982) mark the division between mobile 
and immobile water at eight kPa soil moisture tension. 
Zimmerman et a]. (1965) assign the volumetric ratio of 3:1 
for immobile:mobile water. A theoretical basis for assigning 
immobile and mobile solution pools has not yet been 
presented. 

Solute transport through structured soil cannot be 
adequately described by a longitudinal spreading coefficient 
(Dspr) in the chromatography equation. The numerical 
approach of dividing the soil solution into two phases, or 
peels, anvorder tosaccount for holdback “of solute snethe 
micropore region has been successful in studies of leaching 
in structured soil and the mobile-stagnant pool model is the 
Simplest example of this. 

In summary, the chromatography differential equation 
provides a theoretical framework for solute transport 
Studies in soils, however, the equation must be adapted for 
leaching studies on structured soils. The two solution pool 
concept has proven useful in describing transport under 
these conditions. 

In this study the interactive mobile-stagnant phase 
model of Addiscott (1977) is being tested in simulation and 
field experiments on an Alberta Gray Luvisol. Modifications 
to account for variation in 6 and v with depth and time, and 
for the effects of frozen soil must be incorporated to 
describe the conditions under study. The model is also 


extended to describe transport of adsorbed solutes. 
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II. MATERIALS AND METHODS 


A. SIMULATION MODEL - DEVELOPMENT 

The simulation model described in the following 
sections is based on a numerical solution of the 
chromatographic equation. The soil column is divided into 
discrete layers and water flux between the layers is defined 
by a water balance equation. The solute transport equation 


is then solved uSing rectangular integration. 


1. Moisture Flux 
Moisture flux through the soil column is defined at 


each layer by the following water balance equation: 


WL = WA - storage — 625), 
where WA = water entering layer I from above, 
WL = water displaced from layer I by WA, 
= WA entering layer I+1. 
storage = water absorbed by the layer,I. 


WA added at layer one is determined as the difference 
between the water input and evaporation (WA = W - E). 

This equation is similar to that used by Burns (1974) 
in his layer model of solute transport in non-structured 
soils. The storage available in each layer is defined as the 
difference between field capacity (33 kPa soil moisture 
tension) and the water present in the layer. Water will 


accumulate in each layer until field capacity is reached, 
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after which an amount of water, WL, defined by the water 
balance, will be displaced into the next layer. In terms of 
the symbolism used in the model available storage is 
determined as: 


XSWMCP (I) 


WMCAP(I) + WR(I) - WT(I) (260% 


where WMCAP(I) 


storage capacity of the mobile water pool, 


WR(I) = storage capacity of the stagnant water pool 
WT(I) = total water present at time, J. 
I = layer number. 


Total water present in layer I can be expressed as the 
sum of mobile and stagnant water present in the layer: 

WT(I) = WR(I) + WM(I) R271)5, 
where WM(I) = water present in mobile pool. 

In its present form the model does not deal with solute 
transport at 6 less than the capacity of the stagnant pool. 
The amount of water present in the stagnant pool (WR) is, 
therefore, constant and storage capacity can be expressed in 
terms of WMCAP and WM only: 


XSWMCP(I) = WMCAP(I) + WR(I) - WT(I1) 


WMCAP(I) + WR(I) - (WR(I) + WM(I)) 


WMCAP(I) - WM(I) K269 


The relationship between WA and XSWMCP(I) defines three 
different leaching conditions: 
17) WA < XSWMCP(I) in which case all of WA is 

absorbed by the available storage capacity in layer I; 
ii) WMCAP(I) = WA > XSWMCP(I), in which case an amount 


of water, WL = WA - XSWMCP(I) is displaced into layer 
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I+1; WA is then reset to WL and applied to layer I+1; 
iii) WA > WMCAP(I), in which case all of WM(I) is displaced; 

WM(I) is reset to WMCAP(I); WL = WA - XSWMCAP(I); 

and WA is reset as in ii) and applied to layer I+1. 
Slow leaching is defined by cases i) and ii), while fast 
leaching is defined by 111i). The differences between the two 
will be discussed in the following section. 

The program parameter M defines the size of the soil 
system. Layer M+1 1s regarded as a sink, so that water, WL, 


displaced from layer M is lost to the system. 


Frozen layer routine 

A frozen layer is regarded as an impermeable barrier to 
water movement so that infiltrating water encountering a 
frozen layer must be stored above it. This requires 
definition of a third water pool, WGRAV(I), or gravity 
water, and its storage capacity, WGRCAP(I), and available 
Storage, XSWGCP(I). The relationship among these terms is 
Similar to that defined for the mobile water pool: 

XSWGCP(I) = WGRCAP(I) - WGRAV(I), 
where WGRCAP = total porosity - WMCAP - WR. 
The available storage in each layer, I, above the frozen 
layer is determined as: 

storage = XSWMCP(I) + XSWGCP(I) (2:9. )%. 

The frozen layer number, FLNR, 1s read by the 
Simulation program at each time step, before the water 
inputs are read. If the soil is completely thawed a value of 


0 is entered. If FLNR moves downward between succesSive time 
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steps then the amount of WGRAV above FLNR is redistributed. 
WA is set to WGRAV(FLNRT-1), where FLNRT is the previous 
value of FLNR and FLNRT-1 is the next higher layer. WA is 
then transported downward according to the water balance 
equation (25). The program then moves upward though the soil 
column redistributing WGRAV contained in each successively 
higher layer. 

When the program encounters a layer, I, such that I=0, 
or WGRAV(I) = 0, redistribution ceases and water inputs are 
read for the current time step. If a frozen layer causes 
ponding of water above the surface layer the amount of water 
ponded, WPOND, is added to WA applied to layer one at the 


beginning of the next time step. 


2. Solute Transport (non-reactive solutes) 

Transport of non-reactive solutes through heterogeneous 
media can be written (neglecting spreading effects) in terms 
of the classical chromatographic equation as: 

0(O°c)/dat = - d(v-c)/dz G30). 
Both 8 and v can vary with depth and time. 

A finite difference equation parallel to (30) can be 
written as: 

A(O-ey/At = FRAT hace 7 Az ck be 
The left side of equation (31) can be expanded to: 

A(@c) /At) =) GAc/Atet+ocA8/Ate +o AcA8/At 

If changes in c and 8 with time are small relative toc 


and 6 then the last term can be neglected. The values of the 
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terms outside the difference operators can be set to: 
d>aver= S(6(tyz) et 0 Ch +4°2) )<) and 
C-avyev=tst(c (hy Zz yeteot ttl. 2) ye 

Then: 


A(@-c) = O-aveAc + c-aveAé 


O-avelolemiezy re eclt zat c-aevel0(tez) = <0ttr 1-2) ] 


et tilez ie O Gta lez) es GA tea 0 GE az.) 
Therefore the left side of equation (31) can be written as: 
PECs) Ae wam (ete) SO Gti 7 eeme (tz jdt) 2 i) /ACUS2), 
The right side of (31) can be expanded similarly (if Av 
and Ac are small relative to v and): 
A(v-c)/Az = (v-avedc + c-avedAv)/Az, 
WieeewpoaVer= succes (t,2- 1) )° 
Ceave = sett 2) + cltyzesiod: 
Re sveGo, Zym-ucét , z= 10% 
NP ERR as OFAN) 5 iO Oe wea fen 
Therefore the right side can be written as: 
SMGyrcywAzZG=furhetb, zyrottez)— c(t,271)-v (ty, 2-4) Az, 933)". 
Each of these expressions, (32) and (33), can be 
written in terms of solute masses, M, and water volumes, V: 
c = M/Vw;: 


AV/(A-At); 


VU 


6 Vw/Vb = Vw/(A-Az); 


where, Vw = volume of water in layer of thickness Az, 


and; Vb«=stotal volumesof»plateiof«thickness*Az; 
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Therefore: 


Eon ais Lp Ctaaley Zeina co) eM Uta. VW Ct Hanae) 
At “vw Cah sz AAz Vw(t,z) AAz 
aM yl Zac ieee cir 
AAzAt 
and: 
Rus. S A Gt pai sAvat, ze weiMCenzsi) cAVOE, Mabe 
Az ae z ) AAt Vw(t, z-1) AAt 


= 1 Geeta cell ts a EoM(tyzrisavilt veri) 
AAzAt vw(t,z) Vwcteeoije ae 


Putting L.S. and R.S. together: 


MGCP zo = MCE 2) = (SM (tz) PAV Ct, Zz)! FOM Cty zen Ave ze 14) 
vw(t,z) Vw(t,z-1 
The symbolism used in the simulation model is as 
follows: 
I = layer number (vertical coordinate); 


J = time (coordimate((t/At) 


STRGE) M(z,t) = total dissolved solute mass in layer I; 
WECl jes Vubzathusttotabavatersinelayer lh: 

WM(I) = water in mobile water pool in layer I; 

SM(I) = solute mass in WM(I); 

WR(I) = water in stationery pool in layer I; 


SR(I) solute mass in WR(I); 


WGRAV(I) = water in excess of field capacity, stored above 


a frozen layer; 


SGRAV(I) solute mass in WGRAV(I). 
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WA 


water entering layer I at time J; 


WL 


water dispaced from layer I at time J. 
The rationale behind the two-pool system is explained in the 
literature review. 

WA 1S determined at each layer, I, and time step, Jd, 
from the section of the model describing water movement 
(section ILI. Atts, above). 

Equation (34) can be written in terms of the symbols 
defined above for non-reactive solutes: 

ASTCL) “= WASSTCI— iW) /WT( I= 1) 95>" WLeSt (lL )/Wr (1) (35). 
The solute concentrations in each pool are equilibrated 
within each layer at the completion of each time step so 
that: 

SM(I)/WM(I) = SR(I)/WR(I) = SGRAV(I1)/WGRAV(I) = Sst (1) /wT(t) 
(36). 

Since movement of solute occurs in the mobile phase 
only, equation (35) can be written as: 

AST(I) = WA*SM(I-1)/WM(I-1) - WL*SM(I)/WM(I) (35a). 
where the first term on the right side of the solution 
defines the amount of solute added to layer I from the layer 
above (I-1) and the second term defines the solute lost from 
I to the layer below (I+1). It can be seen from (35a) that 
the mass of solute transported between layers is 
proportional to the volume of water moving between layers 
and the solute concentration in the water. Equation (35a) is 
equivalent to the slow leaching algorithm used by Addiscott 
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For redistribution of SGRAV following movement of a 
frozen layer I is set at FLNRT, WA is set at WGRAV(I-1) and 
(35) is used to determine solute flux between layers. The 
program moves up, decrementing I, until WGRAV is exhausted 


orjthessoit surface tsdreached* 


Fast leaching 

Equation (35a) describes solute flux for the cases in 
which WA < WMCAP(I) (slow leaching). When WA > WMCAP(1) 
(fast leaching), however, the solute added to layer I is not 
given by (35a). In this case the amount of water added to I 
is equal to WMCAP(I) and SM(I) is determined as 
WMCAP(I)*SA/WA, where SA is the solute mass in WA. The water 
added to layer I+1 is determined as WA = WA - XSWMCP(I) and 
solute added to I+1 is determined as SA plus solute 
displaced from I, less solute left in I. This method for 
resetting SA results in mixing at each layer and differs 
from Addiscott's (1977) fast leaching routine which pushes 
solute through as a plug to the layer n, such that 
ZIWMCAP(I) = WA, and then sums the water and solute 
displaced and applies them at layer n. This results ina 
rather discontinuous mixing of infiltrating solution with 
displaced solution. 

In summary the simulation program works as follows: 
by) volume of water entering layer I is determined from 

the water flux routines; 
ii) mass of solute entering layer I (in WA) is determined 


from the appropriate leaching routine; 
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iii) WL is determined from the water balance equation; 

iv) mass of solute leaving I in WL is determined; 

v) ST(1)@iS reset iaccosrdingttoi equdtionwi35) ,f (35a) 4 
or the fast leaching routine; 

vi) solute mass is redistributed between mobile and 
Stationary pools (SM and SR, respectively, and SGRAV, 
if necessary) according to equation (36); 

vii) when ST(I) has been determined for all I the program is 
updated (J=J+1), WA at I=1 is determined as the net 


inéaltrationgatfd=J+1) andsadstohvi) Ware) repeated) 


Evaporation 

The Simulation model adopts the approach of Burns 
(1974) in that it allows the soil to dry out to a moisture 
content defined as the “evaporation limit", 8el. It extends 
the model of Burns to describe transport of volatile 
solutes. 

The net infiltration or net evaporation is determined 
at each time step as WA = W - E, where W is the 
precipitation or water input and E is the potential 
evaporation. When net evaporation occurs the uppermost layer 
can dry to 6el. Any evaporation in excess of this amount 
must be replenished by capillary movement (in the liquid 
phase) from below. Water lost from the surface, WE(1), is 
set equal to net evaporation. If net evaporation exceeds 
WT(1) then WE(1) is set to WT(1) and evaporation 
requirements in excess of WT(1) are stored. The evaporation 


routine is then repeated as many times as necessary to 
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Satisfy the evaporation requirements (-WA). If the amount of 
water (WT(1)) remaining in layer one is less than 6el this 
deficit ((6el - WT(1)) is replenished by water from layer 2, 
so that water lost from layer two is defined as: 

WE(2) = @e1(1) - WT(1). 
The equation defining water loss induced by evaporation for 
i, eoione 1S; 

WE(I)wer dele oe wre ie)s 
This value is assigned until net evaporation loss is 
Satuist ied eithatils Oirormellviyasuchothatrwittiot <del (ie 

The mass of solute lost from each layer, SE(I), is the 
product of the water lost and the solute concentration in 
the solution phase: 

SE(I) = WE(I)#ST(1I)/WT(I). 
Water and solute mass balances are written for each layer 


for whieh Loss toecurs: 


AWT (I) - WE(I) + WE(I+1); 


AST(I) = SGI), ae Soa 


W 


After WT(I) and ST(I) are reset solution within each layer 
is equilibrated according to equation (36). This achieves a 
continuous mixing of solution as it moves upward, similar to 
that which occurs during slow leaching. 

This model extends that of Burns (1974) to the 
description of behavior of volatile solutes, such as ammonia 
and many pesticides. Volatilization losses of solute mass, 


SV, occur from the first layer only and are approximated as: 
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SV = SE(1) = WE(1)*ST(1)/WT(1). 

AS a working approximation the evaporation limit of 
this model is set at the 200 kPa soil moisture tension water 
content. Water in excess of this amount can be lost to 
replenish moisture deficits above. A moisture deficit is 
defined as any 8 less than 6 at 200 kPa. The soil is not 
allowed to dry to tensions greater than 200 kPa. For arid 
solls a lower value for the evaporation limit may be 
neccesary. Alternatively, it might be desirable to allow 6el 
to vary, depending on the overall moisture status of the 
soil profile. 

The evaporation routine was included to make the model 
more general. Due to time constraints no attempt has been 
made to determine a more appropriate value of 6el. The 
definition of SV is valid for tritiated water but may not be 
for other solutes. Further theoretical and experimental 


investigations are necessary to improve this routine. 


3. Transport of Adsorbed solutes 

Partition of solute mass between solution and adsorbed 
phases in soil can be described by a linear isotherm for 
lindane and other pesticides: 

a = Kads-c. 

Changing to the symbols and dimensions used in the 
Simulation model involves the following relationships: 


a =) 07 (DbAZ.)® Fand:, 
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CRE STW); 
where Q = adsorbed solute mass (per layer) per unit 
cross-sectional area; 
Db = dry bulk density; 
Az = layer thickness; 
ST = dissolved solute mass (per layer) per unit 
cross-sectional area; 
WT = total water (per layer) per unit area. 
The adsorption isotherm can be defined in terms of Q 
and Sik: 


Q 


(DbAz-Kads)ST/WT 
= B-ST Se) 
where B is the effective partition coefficient defined by: 

B = (DbAz/WT)Kads. 

The adsorbed solute mass, Q(I), in layer, I, can also 
be expressed in terms of the total solute mass, MS(I), and 
total dissolved solute mass, ST(I), in layer, I; 

OCT) S=aMSOIie-aSh(t) (38)x. 
Therefore: 

SE (i)e= MS(1)7 (A+B) oe 
and, 

Q(I) = B*MS(1I)/(1+B) (40). 

A change in MS(I) resulting from mass flow into and/or 
out of I can be defined by an equation parallel to (35): 

AMSOiijje= SAcsySbhe 
where SA = solute mass entering I 


WA*ST(I-1)/WT(I-1)+3 
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and SL = solute mass displaced from I by WA 
= WL*ST(1)/WT(I). 


Any change in MS(I) causes a change in Q(I) and ST(1) 


so that: 
ST (1 jee MSi(T)e, / (+B) 
=. ¢ MSs) Rees AMS( INO" (aeHB)) (39a); 
Q(I); = MS(I) ;/(1+B) | 
=t BUMS (ie) ano ae AMS) )eAtiHB) (40a), 
where 24 and ; denote successive time steps 


Equations (39a) and (40a) define the partition of 
solute mass between dissolved (ST(1I))° and adsorbed (Q(1)) 
phases at layer, I, for solution moving through I. The 
Simulation program functions as follows: 

(iu) change in solute mass in layer, I, is determined 
asl AMS({i®@ = GSA FeiSL: 
(ii) ST(I); is determined from (39a); 
(iii) Q(1), is determined from (40a); 
(iv) J is updated and continue. 
For non-adsorbed solutes Kads, and therefore B, are set 


abo Oy sor thatwST( ism =eMSKi)echedt AMSs andi Old )O=f De 


4, Pesticide Decay 

Decay of solutes can be treated as a sink term. 
Biological degradation of many pesticides can be treated as 
a first order decay reaction at low concentrations 
(Hiltbold, 1974). Nash and Woolson (1967) determined that 


loss of lindane in soil followed a first order decay rate 
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law. For first order reaction kinetics the sink term can be 
expressed as: 

dS/dt-="-r14s (41), 
where S = the amount or concentration of pesticide, 

r = decay rate constant (1/t). 

Using this sink term involves subtracting a constant 
fraction (r) from the total amount present (MS) at every 
time step. It may be argued that adsorbed solute is 
protected from decompostion and that only the dissolved 
solute 1S Subject to decay. If this is true then the decay 
rate equation would apply to ST (dissolved solute mass) 
rather than MS. The decay rate constant, r, given above is 
an averaged value which defines the overall decay rate in 
all phases and, therefore must by applied to MS rather than 
STambor bindanere rsdlesskthan 0.001/d (Hadit bol dt rikova)e.so 
that 1/1000 X ST is lost each day. 

If time steps are small enough then AS = -rAt:S. The 
analytical solution for (41) is straightforward, however: 

ShitAt) a= iSttiaexpl-—cAtix 
Decay of solutes can be effected in a finite difference 
model by multiplying the solute mass by the factor exp(-rAt) 
at each time step. 

For the period of the experiments performed lindane 
decay was not a factor in solute distribution and a 
decomposition term was not included in the simulation 


program. 
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B. FIELD EXPERIMENTS 


TeSo rh 

Field experiments were performed on untilled barley 
stubble on Breton Loam, an orthic Gray Luvisol. The 
experimental site was located on the University of Alberta 
Breton Plots, two km southeast of the village of Breton, 
Alberta, NE-25-47-4-W5. The soil was formed on medium fine 
textured till and posseses moderate internal drainage 
(Howitt, 1981). 

The physical and chemical properties of this soil were 
described by Howitt (1981). It is characterized by a clay 
ena chedest whorizon 018 — 60cm) with a stronmq blocky 
primary structure and weak prismatic secondary structure. 
The Ap horizon is low in organic carbon (1.4 %) and 
possesses a weak platy structure (Crown and Greenlee, 1978). 
The BC horizon (65 cm -) has a weak coarse blocky structure. 


Sorl characteristics are isummarizedtim itabregii.1. 


2. Solutes 

The solutes used in the field experiments were *HOH 
(tritiated water) and ‘‘C-labelled "lindane". Tritiated 
water was treated as an infinitely soluble, nonreactive 
solute. Lindane (y-hexachlorocyclohexane, also known, 
incorrectly, as y~benzene hexachloride or y-BHC) is an 


insecticidal ingredient in some seed treatments used in 
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Western Canada. It 1S a Sparingly soluble, strongly 
adsorbed, and persistent chemical and, therefore, assessment 
of its mobility in soil is an environmental concern. 


Chemical properties of lindane, pertinent to this study, are 


duven. umitableti®. 2% 


Table II.1. Characteristics of Breton Loam (Orthic Gray 
Luvisol) 


Particle size analysis 


Horizon. Morqs.cl%) %sand 25 %clay 
2AD Ie 33 55 12 
eB 0.4 26 40 So 
BG Ora 30 42 28 


t—Howi1tet ©6198 1,) 
2-Crown and Greenlee (1978). 


Table II.2. Chemical properties of lindane. 


*AGdsorption’ *Kd 


‘aqueous *‘montmor- ‘BretoneL= tenalr life 
IUPAC name Solubriney sconce Ap Bt in soil 
y~hexachloro- 10 ppm aa te 5 2G 


cyclohexane 


'-Green (1974), 
PH lepold. (1974), 
°-Kd expressed as (umoles lindane/kg adsorbent) /(umoles 
lindane/l solution) for equilibrium solution concentration 
of: 10 umoles/l for ‘montmorillonite, and 

4 umoles/l for "Breton Loam. 
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3. Experimental Design 

Tritiated water and solutions of °*HOH and ‘'*C-lindane 
were each applied to the Ap and Bt horizons. The isotope 
solutions were mixed with air-dry, ground and seived (2 mm 
mesh) soil in the laboratory, using a pastry blender. In the 
first treatment 6.25 uCi of *HOH alone were added with 250 
Mee Olmuwater toOmeach Oted.>5 KG. Of Ap horizon soiliand 12.0 kg 
Steept nora zonesos«, 

A second treatment of *HOH-'*C-lindane solution was 
atSOumixedawith so1l., In this case 2.0 4Ci ‘“*C-lindane and 
10.0 uwCi *HOH dissolved in 200 ml of water were mixed with 
SaAcueOre 07 kor OreApeandy 1.5 ko of Bt soa i. 

The soil-water-isotope mixtures were applied to the 
soil on Nov.26 and Dec. 6, 1981, after the soil surface had 
frozen. At 12 siteS approximately 24 cm depth by 30 cm 
diameter of Ap horizon was removed and replaced by the 
Ap-isotope mixtures. Another twelve sites were excavated to 
the upper surface of the Bt horizon and Bt-isotope mixtures 
were added, lightly tamped and covered with Breton Loam Ap 
norizony soil; 

The isotope mixtures were laterally contained within 30 
cm diameter open cylindrical heating ducts pressed into each 
respective horizon to a. depth of two to three cm. The 
heating ducts marked the sites for sampling. The 24 points 
of application of the isotope mixtures were spread out 
evenly over a rectangular plot area eight m X 14 m, as shown 


in figure II1.1. Gamma-densitometer access tubes were 
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= 10 uCi 30H andy Ze uGi See sidindane applied to Bt horizon (17-25 cm). 


+ = Site of itracér application. 

Rea Gi oo. UCL 3508 applied to Ap horizon (0-2 cm). 

B= "6.25: uCt 30H applied to Bt horizon (17-20 cm). 
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x x = Site of gamma densitometry twin probe access tubes. 


Figure II.1l. Field experiments, plot layout. 
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inserted near the N-E and S-W corners of the plot area. 

The plot was sampled between May 5 and May 12, 1982, 
after the snowpack had completely melted. Two 4-cm diameter 
columns were taken from within each cylinder to depths of 
100 to 120 cm using a hydraulic powered core sampler. The 
columns were sectioned into 2.5 cm layers from 0 to 75 cm 
and five cm layers from 75 to 120 cm depth. Solute 
activities were measured by liquid scintillation counting as 


described in section 4.b), below. 


4, Field Measurements 
a. Infiltration/Evaporation 
i) y-Densitometry 

Wet bulk densities, Dt, were determined by the 
y-densitometry method using a TroxlJer Model] 2376 two probe 
density gauge. The method and theory of operation are given 
by van Bavel (1958), Smith et a]. (1967), and Troxler 
Electronic Laboratories (1976). Dry bulk densities, Db, and 
volumetric moisture contents, 6, were also determined from 
data obtained from y-densitometry, as described below. 

Net infiltration or net evaporation were determined as 
the net difference in 8 between successive 86 profiles 
determined by rectangular integration, as illustrated in 


figureml KIM Stinwehbe next-chapter. 


ii) Snow survey 
Total potential infiltration was determined by a snow 


survey on Mar.25, 1982, after the last precipitation event 
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before the field experiment ended. Snow cores were taken 
near the four corners of the plot area usSing a tube of known 
diameter; snow was weighed and height equivalent of water 
determined as the weight over the cross-sectional area of 


the tube. 


b. Volumetric Moisture Contents, 6 

Volumetric moisture contents were calculated as the 
difference between Dt and Db. Ten 8 profiles were measured 
between Mar.9, 1982, before spring thaw began, and May 15, 


1982, after sampling was completed. 


c. Depth of Frost 

Frozen layers were determined by probing with a 
Sampling tool. It was found that frozen layers determined a 
this way corresponded approximately to bulges in the 


difference between successive 6 profiles. 


5. Laboratory Analyses 
a. Soil Physical Properties 
i) Dry bulk density, Db 

Dry bulk densities, Db, were calculated from wet 
densities, Dt, and gravimetric soil moisture contents, w, 
according to the equation: 

Db = Dt/(1tw). 
Gravimetric moisture contents were determined from two 
Samples taken within one m of the access tubes in Nov. 1981, 


and one composite sample taken from between the tubes in 
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Sept. 1982. Samples were removed from the field in sealed 
moisture tins and dried at 105°C. w was calculated as weight 
of water in the samples divided by the weight of the 
oven-dry soil. Bulk density profiles were determined on 

Nov. 20, wide wand Sepr 216, 1962, "and the mean Db calculated 


at 2.5 cm increments to a depth of 110 cm. 


ii) Porosity, a 

Porosities were calculated from Db and particle 
densities, Dp, as: 

am=enia> DbZDp* 
Particle densities were determined by the method described 
by Blake (1965), using volumetric flasks. A porosity profile 
was determined for the same depth and resolution (110 cm and 


2.5 cm, respectively) as for Db profiles. 


iii) Moisture retention curve 

Soil moisture retention curves were determined by the 
pressure plate method described by the manufacturer of the 
apparatus used (Soilmoisture Equipment Corp.). Intact soil 
cores were used in order to preserve the pore structure. 
Steel rings, 3.0 cm high by 4.8 cm in diameter were pressed 
into the soil, dug out and trimmed with a trowel. The same 
cores were used for moisture determinations at all tensions. 
Weights of moist soil were measured between runs and oven 
dry weights measured after the last equilibration. 
Gravimetric moisture contents were determined for Ap (8 


samples), Bt (7 samples), and BC (4 samples) horizons at 5, 
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Sam 1007 e200 fFes00wandeis00 kPacmoisture stensions;bin ithe 
order given. Volumetric moisture contents at these points 
were obtained from the gravimetric according to the 
equation: 


8 = Db-w. 


b. Chemical Analyses 
i)Adsorption 

Adsorption rates and isotherms were determined for 
lindane adsorption to Breton Loam Ap and Bt horizon soil. 
Experimental methods were similar to those of Kay and Elrick 
(1967), except ‘'*C-labelled lindane was used and measurement 
was by liquid scintillation counting (LSC) rather than by 
gas chromatography. 

Two gram soil samples, air dried and ground to pass a 
one mm mesh seive, were Shaken in centrifuge bottles with 50 
minohl' “Gelindanessolutioncof known concentration fand 
activity. Samples were then centrifuged at 27,500 X gravity 
for 30 minutes. ‘'*C-lindane activity in solution was 
determined by LSC of two ml of the Supernatant. Activity per 
unit volume was converted to lindane concentration by 
dividing by the specific activity. The amount of adsorbed 
lindane was calculated as the difference between the 
original and final amounts of dissolved lindane, less 
lindane adsorbed to the centrifuge bottles, determined from 
the isotherm for lindane adsorption to the bottles. 

Adsorption isotherms were determined by shaking soil 


samples in incremental solutions containing 0.005 ppm to 2.5 
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ppm lindane for two hours. Adsorption isotherm experiments 
were duplicated and the adsorption coefficients determined 
as the slope of the regression line of adsorbed versus 
dissolved lindane for the combined sets of data. 

Adsorption aves were determined by shaking triplicate 
Soil samples in 1.5 ppm lindane solutions for time periods 
ranging from 0.25 to four hours. Adsorption equilibrium was 
attained in less than one hour for both Ap and Bt horizon 


Samples. 


Liquid Scintillation Counting (LSC) 

Isotope activities were measured by liquid 
Scintillation counting in Fisher ™Scintivers scintillation 
cocktail on a Searle Isocap/300 System after dark adaptation 
for two hours. Counting efficiencies were determined from 
Sample channels ratios (SCR) as described in the instrument 
manual (Searle Analytic Inc., 1974) and by Vose (1980). Net 
counting rates (sample cpm - background cpm) were converted 
to actual counting rates (dpm) by dividing by the 


efficiency. 


ii) Field sample assay 

Activity of field samples was counted directly in the 
soil samples by a method similar to that described by Lavy 
etealro(1972)). sXanonsgellingiscintillation cocktari (table 
II.3) was used so that soil in the mixture would settle to 
ensure that fluorescence produced by scintillations would be 


transmitted through the cocktail. It was determined in 
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Separate experiments that addition of 15% methanol (by 
weight) to the scintillation cocktail was sufficient to 
dissolve up to two ml water in 16 ml of cocktail at the 
operating temperature of the scintillation counting system. 
The cocktail used could, therefore, dissolve all the soil 
solution in three gram samples of gravimetric moisture 


conten ts 7M w)viups t oF F2e. OP. 


Table 11.3. Preparation of liquid scintillation cocktail for 
direct counting of soil-borne ‘*C-lindane and *HOH. 


Component Amount 
dioxane 1000 ml 
xylene 600 ml 
methanol 400 ml 
W236) 15 gm 
*POPOP 0 4187.52 qm 
naphthalene 240 gm 


'-Scintillants added in a concentrate preparation, Fisher 
™Scintiprep 2, in which PPO and POPOP are dissolved in 
xylene. 

Three grams of moist soil sample were added to 16 ml of 
the scintillation cocktail. Vials containing the mixtures 
were shaken for one to three minutes on a vortex mixer. 
Vials were then shaken for one hour at three cps ona 
reciprocating shaker, allowed to settle in the dark for two 
hours and counted by LSC. 


Sample accountabilities, determined as the net counting 


rate in moist soil standards over the known activity (in 
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dpm) in the standards were determined and an accountability 
venous SCR curve plotted (figure 1Il.2/and figure J1.3)< 

Accountabilities for the direct counting of soil-borne 
‘#C-lindane were determined to be 55% to 61% for Ap samples 
and 68% to 73% for Bt samples. These were on average 5% 
lower than efficiencies determined from the efficiency 
versus SCR curve. The accountability versus SCR curve was 
approximately parallel to the efficiency versus SCR curve 
determined for quenched '‘C standards (figure II.2). 

The direct counting method was extended to meaSurement 
Of activity of “HOH in soil. Aceountabilitiess soid—borne 
7HOH were 15% to 20%, and were consistently 5% lower than 
efficiencies determined from the standard efficiency versus 
SCR curve. The accountability curve for soil-borne *HOH was 
also parallel to wine [totes Qued versus SCR curve for 
soilless standards in the SCR range encountered (figure 
tates 3:). 

Because of the uncertainty involved in extrapolating on 
the accountability curves the standard efficiency versus SCR 
curves were used to determine activities for both 
e4c-lindane and. “HOH. “Samples which contained both %*C and 
7H activities were counted on a dual channel program and 
efficiencies determined from efficiency versus external 
standards ratio (ESR) as described in the instrument manual 


(Searle Analytical Inc., 1974). 
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III. RESULTS AND DISCUSSION 


A. SIMULATION EXPERIMENTS 


1. Relative Sizes of Mobile and Stagnant Pools 

The distribution between mobile and stagnant solution 
pools defines the influence of soil structure on leaching 
patterns in this simulation model. Ratios of WM:WR varying 
between 0:1 and 9:1 were used in a Simulation experiment to 
examine senSitivity to this ratio. The amount of water in 
the mobile pool, WM, was set at WMCAP, the storage capacity 
of this pool, and WT, Az, and water flux density were held 
eonstans(tablesiit.i). (Timervstep, Ati«ftor™thismand all 


other simulations is one day). 


Table III.1 Data used in Simulation experiment ia: Influence 
of size of mobile and stagnant moisture pools (WM and WR, 
respectively) on solute distribution. 


cm/layer 


WATER FLU& 
Az TIME DENSITY 
SIMULATION cm days ml/cm?-d WT WM WR WR: WM 
A 2.54 40 0.5 fea) UB) 0.0 Osa 
B 2.54 40 065 ve.0 Ca Or5 ee 
C 2.54 40 O25 ea) OR Zed Oe a 
D 2.54 40 OFS venO O's! we, act 


The simulation approaches the analytical solution at 


low WR:WM ratios (figure III.1). For WR:WM of 1:1 or less, 


61 


10565. 
. 7 
i G 
4 és 
a 
a 
wv 
. 
© 4 ‘ 
Si 30 


slood tnsapet2 bas elideM 20 4 


apesj2a Sos alidom mkewiad nokF 


ef 


syutzuta2 Lice Ye esoneuglini at onits 


ijsiunte «6 ni Ssaevs S30 te: -2 6a6 i 


neb sul? zetew Sue ,.S6- , TW bine <toog 4 is 
Taz 44. ,Ge3e smgt) fit sides) 3 


an 


foog etesetom Inane 


> 


So 
—. | 


soiiteh. .febom noizvelumia elds ob 


4 


wee 
na atT .St3a3 eid? o2 ¢#ivistoe 


TAI 46 80 Sev . Me feog: i 


(veh 2a0 af snotgatumte % 


beau ‘atad ft. 11% 
bas slides fo si 
citudixze2ib etutlog ac (ylevat 


rrsqgxe noite lumi " a 


5 
ai 
+ 


4 


“ms 


9 2.0 or 2202 = 
} 2.0 Os 2-85 
,! 2.0 of ee. 
Oo €,0 0S $2.8. 


J | : 
ityvisns sat eestizveci9qgs noite Lum: 


‘e 


MWiSW 120% ff, 129 ould} pega 


4 


62 


Depth (cm) 


—— — —— Simulation A (WR:WM 
Simulation B (WR:WM 


Simulation C (WR:WM 


— — — Simulation D (WR:WM 


Analytical solution. 


0 | .05 . 10 15 
Mass / Total mass added 


Ripire Lily... Simulation experiment la: Influence of relative sizes 
of stagnant (WR) and mobile (WM) water pools on distribution of solute 
added to the surface layer (0 - 2.5 cm) of a structured soil after 


simulated leaching by 20 cm of water. 
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Simulated solute distributions are identical and the 
numerical solution is approximately equal to the analytical 
solution. Viewed qualitatitively this means that low WR/WM 
values simulate the situation in which structure plays a 
decreasing role in affecting leaching patterns. WR/WM 
approaches zero for structureless media. This is illustrated 
in figure III.1 by the special case of WR:WM of 0:1 
(simulation A). High WR/WM values, representative of highly 
aggregated soils, are illustrated in simulations C and D. 

Increasing the influence of soil structure (increasing 
WR/WM) iS accompanied by both a shift in the solute 
concentration profile toward the soil surface (positive 
skewness) and greater spreading. These effects are produced 
by a fast leaching routine in the simulation program which 
causes channeling of soil solution through the macropores 
when the incremental infiltration inputs exceed the mobile 
water capacity (WMCAP). 

This fast leaching routine reduces the degree of 
equilibration between mobile and stagnant solution pools. 
Where the rapidly flowing solution moves through a region in 
which stagnant pool concentrations are higher than in the 
infiltrating scolunion, diffusional mass transfer from the 
stagnant pool is reduced and less solute mass is transported 
than for slowly infiltrating waters. In regions where 
concentration in the stagnant pool is lower than in the 
moving solution reduced mass transfer between pools results 


in transport of greater mass to greater depths than would 
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occur under low flow rates. 

Thus sthepnet ieffectaciethe fast leaching routinevis to 
cause deeper penetration of a small amount of solute but 
holdback of a greater proportion of the mass of the solute 
at the same time as indicated by the upwardly skewed 
concentration versus depth profile and the the tailing of 
the profile at greater depth. This phenomenon was observed 
by Dekkers and Barbera (1977) in column leaching experiments 
with initially dry aggregated soil. 

The analytical solution in figure III.1 provides a 
check on the accuracy of the numerical techniques used in 
the simulation model. The concentration peaks are at the 
Same position in the soil column for the analytical solution 
and the structureless system. Slightly more spreading occurs 
in the simulation, however. 

The dispersion coefficient, E, used in the analytical 
solution (Kirkham and Powers, 1972): 

C/Cor = Oso Cert (zi teo-r bb) ev (Bt) le eert (2 -vtyy cv (ste C4 i 


Where Zsa 2 = 2.545cm> 


ZG al eka ee Ken |S 


p fluid velocity 


PEAT elon Waren 


E AY ds 


was determined as the numerical spreading coefficient using 
the approximation formula given by Reiniger and Frissel 


(1974): 
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DSprP= b-Az7z2 (42). 
Using this approximate value of E introduces error into the 
analytical solution. Values for E determined from (41) and 
the maximum concentration from the profiles obtained from 
Simulations A or B yield curves that effectively coincide 


with the analytical solution. 


Determination of WR and WM 

Soil aggregation influences leaching patterns because 
of the heterogeneous pore structure it produces. Therefore 
it 1S important to divide mobile and stagnant pools on the 
basis of pore size distribution. The moisture retention 
curve (6 versus soil moisture tension) is one way of 
quantifying pore size distribution, since soil moisture 
tension is a function of pore radius. 

To assign values for WR and WM a moisture tension 
characteristic of the pore size range of each pool is chosen 
and volumetric values assigned on the basis of the moisture 
retention curve. In this way a standard theoretical 
criterion for dividing soil moisture pools can be used. The 
moisture tension used to assign WR/WM should be the same for 
ald SoilshtAddisecttyoil977) .. Addiscottréi977) Ghastusedo200 
kPa tension as a working approximation. 

Sensitivity to the moisture tension value used in the 
allocation of WR:WM was examined using the data in table 
III.2. The upper limit of WMCAP was defined as field 
capacity (@ at 33 kPa soil moisture tension) and WM was set 


at WMCAP. The upper limit of WR was set at 500 kPa, 200 kPa, 
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and 100 kPa for separate simulation runs. Values of 8 were 
determined from a soil moisture retention curve for 
undisturbed soil cores taken from the field. A complete set 
of soil physical properties were determined for this 
experiment by field and laboratory measurements (table 


5 Gd ORNS) ee 


Table III.2 Data used in simulation experiment 1b: Influence 
of soil tension value used to assign WR:WM on solute 
distribution. 


cm/layer 


WATER 

FLUX 
SIMU- Az TIME DENSITY WR: WM 
LATION cm days, ml/cm?-d WT WM WR set at 
B 2.5240 ORS varieS varies varies 500kPa 
C 225m S40 Oats) varies varies varies 200kPa 
D 2.54 40 Orn varies varies varies 100kPa 


The results of this simulation experiment (figure 
III.2) were in qualitative agreement with the previous 
experiment. Choosing a low soil tension to assign WR and WM 
resulted in a high WR:WM ratio, cauSing a greater positive 
skewness in the solute concentration profile (simulations C 
and D). The skewed profiles indicate greater penetration of 
solute and lower net movement of solute mass, as described 
in the previous Simulation experiment. Increasing the 
tension value (simulation B) moves the solute concentration 
toward the analytical solution (A). 

The abrupt change in slope of the solute concentration 


profile at 15 cm depth coincides with the boundary between 
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Figure III.2. Simulation experiment 1b: Influence of soil moisture 
tension values used to assign sizes of stagnant (WR) and mobile (WM) 
water pools on distribution of solute added to the surface layer 

(O - 2.5 em) of a structured soil after simulated leaching by 20 cm 


of water. 
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Ap and Bt horizons and reflects the differences in moisture 
retention characteristics between these horizons. The higher 
clay content of the Bt horizon and accompanying larger 
micropore volume result in a higher WR:WM ratio in the Bt. 
Consequently a greater proportion of solute 1s retained in 
the’ stagnantypeol ine@the Btzthan in thefApOhérizon. This 
feature causes a sharp increase in solute mass held in the 
stagnant pool at about 15 cm depth. 

Such identifiable features can be used to assess the 
influence of structure on field experiments. Existence of 
Similar sharp changes in slope of the solute concentration 
profile observed in field experiments would confirm the 
relative importance of structure in leaching and aid in 
determining the appropriate distribution of WR and WM. Other 
diagnostic features of the leaching system which show up in 
the simulation runs are the position and magnitude of the 
solute concentration peak, and the depth of solute 


penetration. 


2. Rate and Pattern of Infiltration 
High infiltration rates, Simulated using the data in 
table III.3, have an influence on the system similar to that 


of low WR:WM ratios. 
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Table III.3 Data used in simulation experiment nr. 2a: 
Influence of infiltration rate on solute distribution. 


cm/layer 


WATER FLUX 
Az TIME DENSITY 
SIMULATION cm days= mil/cm-—d WT WM WR WA /WMCAP 
A 2.54 80 O25 en he 5 Ome 
B Zoe 40 0.50 ets eon) HO Omer tO) 
C 2.54 20 16 Tee ge) Soe ar eee) 
D 252 8 Pees eee Ove ae Ore 


The quantitative relationship between the size of water 
inputs, WA (WA = flux X At), and the mobile water capacity, 
WMCAP, determines whether solution is channelled through 
soil macropores by the fast leaching routine (WM is set at 
WMCAP). For WA/WMCAP greater than 1, solution moves through 
the soil column with limited equilibration between SR and 
SM. This mechanism J apimarees channelling of solution (figure 
eT UNS Oe 

Increasing WA/WMCAP to two and five (simulations Ceand 
D, respectively) causes an increasing upward skewness in the 
solute concentration profile. The simulation mechanism 
causing this skewness is discussed in the previous 
simulation experiment. Increasing WA/WMCAP has an effect 
Similar to increasing WR/WM. 

Simulation B is the special case in which incremental 
infiltration inputs are equal to WM (and WM = WMCAP) at each 
layer. The effects of this situation on numerical dispersion 
were described by Leistra (1973). Mobile solution pools are 
displaced by exactly one layer thickness for each time 


increment and, therefore, no mixing occurs in the mobile 
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Depth (cm) 


/ / 
/ —_— — Simulation A (.25 cm water/day). 
t 
/ ee THOT eleva Uy POO) ain water/day). 
100 # | 
! _ — Simulation C (1.0 cm water/day). 


eee ee Stim at on D.( 2. Seem wat ese dan le 


0 205 LO ais. 


Mass / Total mass added 


Figure III.3. Simulation experiment 2a: Influence of infiltration 
rates on distribution of solute applied to the surface layer (0 - 2.5 


cm) of a structured soil after simulated leaching by 20 cm of water. 
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pool. The spreading observed in this case is entirely due to 
holdback of solution in the stagnant pool. This is why B 
exhibits less spreading than A. 

For Simulations A and B, WA/WMCAP is less than or equal 
to one, so that slow leaching, which permits more complete 
equilibration between SR and SM, is the only transport 
mechanism. A and B are symmetrical about a peak 
concentration located at a depth equal to the product of the 
fluid velocity and time. These plots approach the analytical 
solution. This implies that for WA/WMCAP less than one the 
System approaches an ideal chromatographic system: that is, 
soll structure has little influence on the leaching pattern 
which develops. This fact can be illustrated with box 
diagrams representing simulation of a simplified system 
(figure III.4). Diagrams A, i) to iv), represent a model of 
leaching of layered soil column. This is a single solution 
HOG model swith vasunliorm, poressi ze. eDiagrams By ih) .to vin), 
represent a layer model of a structured soil possessing two 
solution pools, equal in volume and representing two 
distinct pore size ranges. 

At time, to, 100% of the solute present in both models 
is located in the top layer. In the non-structured system, 
A, the entire solution is mobile. In the structured system, 
B, 50% of the solute is contained in the mobile pool, SM, 
and? o0% tS 91n eEhe stagnant pool, SR. 

If a water volume equal to 25% of the layer pore volume 


moves into layer one then 25% of the solute in layer one is 
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displaced into the next lower layer. If the same volume of 
water moves into the structured system 25% of total solute, 
or 50% of SM, is displaced. The resulting solute mass 
distribution for each system, at t,, given in A ii) and B 
iii), are equal. Repeating the above steps for two more time 
intervals result in the solute distributions shown in A iii) 
and iv) and B v) and vii). Solute distribution is the same 
for each system at each time step. Solute distributions for 
case A can be determined from the analytical solution. 
Distributions for the structured system, B, can also be 
predicted by the analytical solution of system A, provided 
WA is always less than WMCAP. 

When WA exceeds WMCAP the amount of solute displaced at 
any time step in the portent system can never be more 
than SM, whereas all the solute in the non-structured system 
is subject to displacement. The structured system causes 
holdback of solute in the stagnant pool and at the same time 
reduces the effective conducting porosity so that greater 
solution velocity results from the same flux density. Both 
these features can be seen by comparing simulations A and D 
of figure III.3, which represent WA/WMCAP values of 0.5 and 
5, respectively. 

This simulation experiment illustrates the relative 
leaching efficiency of different infiltration rates. Heavy 
rainfall events would transport proportionally less 
Soil-borne solute than would gentle precipitation. Greater 


mass of rain-borne solutes, such as sulfate and radioactive 
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fallout, would be carried to a greater depth by more 
intenSive precipitation events. 

The results of this simulation experiment also point 
out the influence of the time step on simulated leaching 
patterns. Activation of the fast leaching routine depends on 
the infiltration rates. The time step used must be small 
enough to detect short events of high intensity. If large 
time intervals are chosen the information on rates of 
infiltration will be hidden in an infiltration rate averaged 
over the entire time interval. 

Simulation experiments were performed to assess the 
effect of pattern of infiltration on a structured soil 


system using data summarized in table II1.4. 


Table III.4 Data used in simulation experiment 2b: Influence 
of infiltration patterns on solute distribution. 


cm/layer 
WATER FLUX 


Az TIME DENSITY 
SIMULATION cm days WT WM WR miycn- -d 
2.54 40 F O25 Oa ace! Sita ell. 5 


Zod 40 Ore OS Ore5 


130 

2.04 100 1.0 02:25 On 75 ON 
ie 0 

Peake 16 0 O25 OFS Mee 


OQW }Y 


An infiltration pattern similar to figure II1.5 was 
observed in experimental plots by y-densitometer scanning 
during spring melt of the snowpack. To assign variable daily 
infiltration rates a smooth curve was plotted through the 


rates measured; the area under the curve was taken to be 
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20.0 cm of water; and daily infiltration was determined by 
interpolation. 

The Simulation for varying infiltration rates (figure 
III.6, Simulation A) is intermediate between the average 
rate simulated (C) and the maximum rate simulated (D). The 
position and magnitude of the solute peak concentration for 
A bear closest resemblence to D (maximum infiltration rate). 
The poorest fit occurs between A and B (the minimum 
infiltration rate). 

Results of this simulation experiment indicate that 
relatively few infiltration events of high intensity may 
dominate over more numerous events of low intensity as shown 
by the the relative agreement of simulation A to B and D. It 
also indicates that the model is sensitive to the degree of 
resolution of the infiltration data as shown by the 
differences between plot A and plot C, the simulation of 


averaged infiltration rates. 


3. Layer thickness 

The influence of layer thickness, Az, on numerical 
chromatographic models was discussed in Chapter I, under 
"plate theories". Leistra (1973), and Frissel and Reiniger 
(1974) have explained how layer thickness affects spreading 
of the solute concentration band in numerical models of 
solute transport based on chromatographic theory. 

The influence of layer thickness on simulated leaching 


patterns.<(figure:J11.7)i «was assessed for three, different 
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Figure III.6. Simulation experiment 2b: Influence of infiltration 
pattern on distribution of solute added to the surface layer of a 


structured soil after simulated leaching by 20 cm of water. 
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thicknesses using the data in table III.5. 


Table III.5 Data used in simulation experiment 3: Influence 
of layer thickness on solute distribution. 


a em/ layer ot 


WATER FLUX 
Az TIME DENS Ty. 
SIMULATION cm days mi/emn*-d WT WM WR 
A Tees 20 Tag OS Uae 0525 
B Oe 20 0 et O25 Wate 
S 10.16 20 a0 4.0 eal) 220 


For this simulation WM was set equal to WMCAP. 
Simulations A, B, and C exhibit slightly different amounts 
of spreading and net solute transport. The smallest layer 
thickness, 1.27 cm, caused more spreading than did 2.54 cm. 
This seems to contradict the prediction of the analytical 
solution. This anomaly can be explained with reference to 
the previous subsections. Fast leaching in this model 
depends on the relationship between WMCAP and WA. WMCAP is 
directly proportional to layer thickness. A small layer 
thickness will reduce WMCAP and thus increase the frequency 
of use of the fast leaching routine. This results in greater 
spreading of the solute profile and greater positive 
skewness of the solute concentration profile. Slightly 
greater positive skewness can be observed in plot A (layer 
ERIVCKNeESS = 1.27) CMm)sitnan OCCULS- IN} plots: Band sc,. for 
greater thicknesses. 

Differences in plots A, B, andvC resulting from 


differences in layer thickness, are small. Reducing layer 
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Figure III.7. Simulation experiment 3: Influence of layer thickness 
on distribution of solute added to the surface layer ( 0 - 2.5 cm) 


after simulated leaching by 20 cm of water. 
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toacknessnby 1/79, Pirom iid: téecmiato: 142 bremenresults tanta 
Ghange in peak "concentration: of vabouty 10%, anda shift (in 
depth of the peak concentration of about one cm. Reference 
to figure III.1 shows that a change in WR/WM of a similar 
scale (plots III.1.B and III.1.D) reduces peak height by 2/3 
and shifts peak concentration position by 20 cm. Inspection 
of figure I11.4 shows that the infiltration rate also has a 
much greater influence on simulated leaching patterns than 
does layer thickness. 

This relative insensitivity to layer thickness occurs 
because changing layer thickness brings opposite forces to 
bear on spreading. Increasing layer thickness increases 
numerical spreading by increasing the mixing height of the 
system, but at the same time reduces spreading resulting 
from channelling through macropores, affected by fast 
leaching. 

In theory, layer thickness should be selected on the 
basis of a physical experiment. The simulation experiment 
described above indicates that the model is not highly 
sensitive to layer thickness between 1.27 cm and 10.16 cm 


and a guesstimate value may be adequate. 


4, Frozen Soil 
A frozen layer can influence solute transfer by 
restricting infiltration rates. Four -simulations, were run to 


assess the effects of infiltration into a partially frozen 
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soil (figure III.8) using the data is presented in table 


Pole. oO. 


Table III.6. Data used in simulation experiment 4: Influence 
of frozen soil on solute distribution. (Infiltration rates 
are 0.5 cm/day in all cases. Flow rates vary because of the 
influence of frozen layer). 


cm/layer 


Az TIME WM & RATE OF THAW 
SIMULATION cm days WGRCAP WMCAP WR cm/d 
A Zn OH ELS (OAs) OrZ5 DER OURO 
B Zoe 25 O45) One 238) OR (ay 9 TPA 
& 2 Oa kao 0225 C25 Ove Seem Or 
D Dik ame Ove OR ZAS) Wy 783) a 


In simulation A, a perched saturated zone forms above 
the impermeable frozen layer and extends upward above the 
soil surface. Downward movement behind the receding frozen 
layer is entirely due to movement of gravity water (WGRAV). 
This water represents a smaller fraction of the total than 
in the thawed soil because there is more water present in 
this saturated zone. Though equal volumes move at each time 
step in each case, a smaller fraction of the total solution 
present moves above the frozen zone. Therefore a smaller 
fraction of the solute is transported in the partially 
frozen soil. 

For simulations B and C the ice recedes much faster 
than in A. The saturated zones in these cases do not extend 
as far as in A. Therefore solution moving in the unsaturated 
upper regions of systems B and C represent a greater 


proportion of the total solution than in A, thus 
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—— — — — Simulation A (frozen layer recedes 
atjoecn /) day). 


— — — Simulation B (7.5 cm / day). 
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Simulation D (completely thawed soil). 
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Simulation experiment 4: Influence of a frozen layer 


on distribution of solute added to the surface layer (0 - 2.5 cm) of 


a structured soil after simulated leaching by 12.5 cm of water. 
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transporting a proportionately larger mass of solute. 

The other major distinguishing feature of these plots 
is the symmetry in plot A and the skewness in B, C, and D. 
Solution flow rate is reduced by the frost layer in case A 
so that greater equilibration between SR and SM occurs. In 
cases B through D the influence of the frozen layers on flux 
rates 1s reduced so that fast leaching is more dominant, 
producing the greater skewness evident. 

When a frozen layer has the effect of reducing the 
solution flow rates and therby causing greater equilibration 
between mobile and stagnant solution pools, it can cause 
leaching of a greater solute mass than would occur in thawed 
soil with the same infiltration rate. An illustration of 
this phenomenon is discussed in the final section of this 


chapter (figuresili 720): 


5. Adsorption 

Adsorption affects retention of solute in opposition to 
the influences of mass flow and diffusion. Quantitative 
effects of adsorption on solute distribution were simulated 
uSing data in table III.7. Solute was added to the first 
layer of soil (2.54 cm thickness) and equilibrated between 
solution phase pools and adsorbed phase according to a 


linear adsorption isotherm. 
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Table III.7. Data used in simulation experiment 5a: 
Influence of adsorption during leaching on solute 
distribution. (Layer thickness and time step were 2.54 cm 
and one day, respectively, for all simulations.) 


cm/layer 
WATER FLUX 

TIME DENSITY Db Kads 

SIMULATION days ml7em--d WM WR g/cm? ml/g 
A 80 @525 O25 0275 1185 0.0 

B 80 OmZ5 O25 OS yee iene) 

C 80 Go25 QO.25 Gen 5 SU 

D 80 Ure O25 O75 eke: 2 nw 


It 1S assumed that adsorption of the hypothetical 
solute is defined by a linear isotherm and that equilibrium 
between adsorbed and dissolved phases is attained at each 
layer between time steps. Therefore the relationship between 
mass (q) of adsorbed solute per unit mass of soil and 
concentration of dissolved solute (c) can be defined as: 

gq = Kads°c. 

The values of Kads used cover a range of conditions 
encountered in soil. Adsorption of some anions to some soils 
Can be ‘detined by an kads of 0) Kads ot 17,55, ands25 
represent weak, moderate and strong adsorbate-soil 
affinities, as, for example, lindane adsorbed to sand B 
horizon, loam B horizon, and loam Ap (0.M. = 4%) horizon 
soils, respectively. 

The most striking feature of the plots of these 
Simulations (figure III.9) is the low mobility of adsorbed 
solutes. (Solute distributions in figure II1.9 are expressed 
as mass per unit bulk volume rather than concentrations 


because part of the solute present is not in solution 
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phase.) Nearly 100% of the non-adsorbed solute has moved 
beyond the furthest extent of even the weakly adsorbed 
solute. Peak concentration of the non-adsorbed solute 
(Simulation A)”~moved 50 em, while for the adsorbed solutes 
peak concentrations moved ten cm (Kads = 1 ml/g), one cm 
(Kads = 5 ml/g) and zero cm (Kads = 25 ml/g). 

Examining simulation B in detail reveals the reason for 
the immobility of adsorbed solutes. At each time interval in 
the simulation, solute is partitioned between solid and 
liquid on the basis of Kads and the relative masses of solid 
and liquid in each layer. Assuming a linear isotherm and 


Hocaleequilibriumy;  thrs partition 1s “defined iby (39): 


ST = B*MS/(1+B), 
where B = Db-Kads/8, 
Dbi="4dry bulk density sof soil (g/em*): and; 
86 = volumetric moisture content(ml/cm?). 


The solute mass, SM, subject to displacement is defined 
by equation (36) so that: 

SM = (WM/WT)ST = (WM/WT)B*MS/(1+B). 
BoteKads#=51.0 mi/g,, andupb and 0 detined swnotableslil./,. 3B 
is 3.8, while WM/WT is 1/4, so that: 


SM 


ite 


MS/20. 

This means that the mobile solution, WM, transports 
approximately 1/20 of the total solute mass as it passes 
through layer 1, where the solute was originally placed. As 
the solution passes through the next layer approximately .95 


of the dissolved solute mass is retained. This partition 
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Figure III.9. Simulation experiment 5a: Influence of adsorption on 
distribution of solute added to the surface layer (0 - 2.5 cm) of a 


structured soil after simulation leaching by 40 cm of water. 
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process is repeated for lower layers at each successive time 
step so that when the solution passes out of the sixth layer 
(Setamet steps Laterdlitecontainseonkygtd/20)ey ortkéss®than 

3X10°-7 of the solute mass which it carried out of layer 1. 

The above numerical example illustrates that 
concentration of leaching solution is an exponential 
function of the negative product of B and time. This 
explains why the plot of mass vs depth drops off so sharply 
ROrP abGimulath OneDrawithekKads Of 25 ml/g (B = 93.75). This 
analysis illustrates that adsorption is a dominant control 
on leaching of solutes through soil. 

Water flux was earlier shown to affect leaching 
patterns of non-adsorbed solute. Transport of a weakly 
adsorbed solute (Kads = 1.0 ml/g) was simulated at two 
solution flux densities (table III.8) to assess the 


influence of flow rate on distribution of adsorbed solute. 


Table III.8. Data used in simulation experiment 5b: 
Influence of solution flow rates on distribution of an 
adsorbed solute (Layer thickness and time increment were 
2.54 cm and one day, respectively, for all simulations). 


cm/layer 
WATER FLUX 
TIME DENSITY Db Kads 
SIMULATION days nic -d WM WR G7eme ml/g 
A 100 O25 R25 en floss Ae 
B 50 0.5 Oa 5 Ond 1S a0 


The results (figure III.10) exhibit the same 


qualitative characteristics as shown in figure II1.3 
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(influence of solution flow rates on leaching of 
non-adsorbed solutes), that is, the concentration profile is 
skewed toward the soil surface for higher flow rates (B). In 
general the phenomena described in the previous simulation 
experiments performed on non-adsorbed solutes also occur 
during leaching of adsorbed solutes but the relative 


importance of their influence is reduced. 


B. FIELD LEACHING EXPERIMENTS 

Radioactive labelled solutes were applied to the Ap and 
Bt horizons on Nov. 26 and Dec. 6, 1981, as described in 
Chapter II (Methods and Materials). Sites were sampled on 
Maye7 ‘to ilz, (982, approximately 10 to 15 daysvafiter the 
snow pack had completely melted. During the first three days 
of sampling a frozen layer was encountered at 75 to 85 cm. 
Sampling on May 11 and 12 did not encounter frost. Samples 
were segmented into 2.5 cm layers from 0 to 75 cm depth and 
five cm layers: frome75 to 100 cm depth. 

(Distribution of tracer for experiments discussed in 
this section is reported as activity per layer over the 
total activity added, (A/A,.) and is also tabulated in 


Appendix A.) 


Experiment 1: Leaching of Non-reactive Solute (*HOH) 
a. °HOH applied to Ap horizon 
Only one site at which *HOH was applied to the upper 


layer of the Ap horizon yielded measureable activity 
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Figure III.10. Simulation experiment 5b: Influence of infiltration 
rates on distribution of adsorbed solute (Kads = 1.0 ml/g) added to 
the surface layer (O- 2.5 cm) of a structured soil after simulated 


leaching by 20 cm of water. 
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(greater than 100 total counts over background) when sampled 
ieene Spring (figure sll 11). Tritium could have beens lost 
by evaporation and subsurface lateral flow through the soil. 
Evaporation occurred early in the winter and late in spring. 
The plot area was bare of snow for short periods at both 
these times. The site which retained measureable tritium 
activities was covered in a snow drift until after Apr. 26, 
1982 and so was less subject to evaporative losses than the 
other sites. 

Tthe plots were located on a Slope’ of 12% to 5% (Crown 
and Greenlee, 1978) and this may have caused subsurface 
lateral flow, especially when infiltrating water encountered 
an impermeable frozen layer present during most of the 
period of this experiment. (Heating duct which marked the 
location of tracer application did not extend below the 
level of tracer application and, therefore, did not restrict 
subsurface lateral flow.) 

Tritiated water applied to the soil surface at this 
Site penetrated to a depth of 25 cm. Two tritium peak 
activities occur in figure III.11: one at the surface layer 
(0 to 2.5'cm) and one at ten to 12.5 cm. The primary peak at 
the surface may reflect assimilation of tritium by soil 
biomass during several days that the tracer-soil mixtures 
were stored at room temperature before they were applied in 
the field. Tritium in *HOH readily exchanges with hydrogen 
in biological tissue and assimilated tritium would be 


immobilized against both leaching and evaporation. 
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Figure III.11. Tritium activity distribution observed in May, 1982, 
after 3408 was applied to the Ap horizon (0 - 2.5 cm) in Dec., 1981, 


for the one site where tritium activities were detected. 
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AsSimilation of less than 1% of total added tritium 
aCtIVITY whore pCiy ag Soll. woulo account Tor ally or the 


activity detected in the surface layer. 


b. *HOH added to Bt horizon 

Trivinm activity castributron for HOH applied: to -the 
top of the Bt horizon 1S Summarized as the profile of the 
mean activities at each layer sampled (figure III.12). The 
Bt horizon was encountered at different depths and so the 
7HOH was applied at different depths. Depths for this plot 
are therefore expressed as distances from the original 
placement (0 cm). Distances above this level are expressed 
as negative values, distances below as positive values. 

Data shown are the means of ten sites, except at 0 cm 
and 2.5 cm where presence of high '*C activities in the same 
Samples made determination of tritium activity impossible. 
Tritium activities reported for these two levels are the 
means of four samples. 

The peak of the mean *H activites vs depth plot occurs 
at 7.5 cm below the level of application. Measureable 
tratium aCervity owas found at 28 cm below, tandvat .15.cm 
above, the level of application. The lower part of the curve 
shows the typical bell shape also seen in the simulations 
presented in section A. The amount of tritium activity 
present in the upper region of the curve did not appear in 
any Simulations. Upward movement of solute could be a result 


of evaporation or diffusion. 
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Figure III.12. Mean tritium activity distribution observed in May, 
1982, after 3408 was applied to the upper surface of the Bt horizon 
in Dec. 1981. Depth coordinate given as distance above (negative 
values) or below (positive values) level of application. Standard 


deviations shown. 
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Measurement of the soil moisture profile by 
y-densitometry detected 2.7 cm of water loss by evaporation 
between Apr. 26 and May 15, after sampling was completed. 
Successive @ profiles on Apr. 29 and May 3 showed water 
depletion to a depth of 45 cm (figure III.13). The upward 
displacement of *HOH indicated by figure III.12 could have 
been caused by water being drawn up through capillaries in 
response to a soil moisture potential gradient induced by 
evaporation from the soil surface. 

Extra diffusional mixing may have occurred in a zone of 
high moisture content perched above a frozen layer. 
Measurement of the @ profile detected a zone of high 
moisture content which showed up as a bulge in the 
difference between successive 8 profiles. Physical probing 
of the soil showed that this bulge often coincided with the 
depth of thawed soil. Gillies (Gray et al., 1970) observed a 
distinct upper saturated zone of thawed soil above a frozen 
zone during infiltration of snow meltwater. 

A frozen layer was present for nearly the entire period 
of the experiments. Diffusion of solute in the zone of high 
moisture content above this layer would be mainly upward 
because of the ice barrier below. Shearer et a]. (1973) 
stated that diffusion rates were approximately proportional 
to 64264, .This.celationship implies that diffusion would be 


greatly enhanced in this zone. 
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Experiment 2: Leaching of Adsorbed Solute ('‘C-lindane) 

Adsorption isotherms were determined for lindane 
adsorbed to Breton loam Ap and Bt horizon soil (figure 
III.14). Points plotted show the results of two separate 
experiments in both cases. Linear isotherms describe the 
data with a high degree of correlation, as determined by 
regression analysis (r? = 0.995 for Ap isotherm; r? = 0.92 
for Bt isotherm). A slight downward concavity is evident in 
the isotherms and fitting the data to a non-linear 
Freundlich 1osthesm (q =. Kads-c'”’"™) “improved the £16 
Mero nally wr ame 3995 hor Ape lsotherm: r2 (=s0 997eronest 
isotherm). 

The Kads for lindane adsorbed to Breton Loam Ap and Bt 
horizon samples were determined to be 12.6 and 5.1 ml/g, 
respectively, from linear isotherms. These values represent 
Strong and moderate adsorption interactions. Kay and Elrick 
(1967) determined Kads for lindane adsorbed to soil with 
higher organic content than Breton loam to be between 17 and 


23 ml/g. 


a. '*C-lindane applied to Ap horizon 

For lindane applied to the surface of the Ap horizon 
98.0% of measured activity was found in the zone of 
application (figure III.15). The observed distribution shows 
low mobility of adsorbed solute and serves to emphasize the 
influence of adsorption on solute leaching. The data 
presented bears strong resemblance to the simulation 


represented by figure III.9D, for movement of strongly 


ae): 

(epsbaif-3'*) stufoe | 
snsbaii a6} 8s ae 
siupii) [fde2 nosiod gat 
etaieqgs ovd 16 2tluasz od woe Salo 
od? sdixszesb. emisryoat gsaentd 29265 ae 
vd banimis245 26 ,sefdsiayies Io ss19sh 
O = *y ymssditeet GA 363. 2e.Ore *F) efay! 
inshive 2i gtiveonés bisenwob siptia A + (nm 
“t6enili-noa 6, 62 sieh ea? parz2tt bas: 

ti? 9d3 bavorgm? ("°"a-ebew = p) GepeRE 


Cer.U = > mies Jon! GA 363 ¢2e8 .0 = ty) 


bacttag 4 


7.4 d. od S9nimeei36e5 sis¥ esiqmae 10S; 


2eulsv o@alT .emtettoei aseail mow? (ylevioee 
; _ 


‘ aa 
m bas prio’ 


a? ~ehOoraves?7i IO Das S20G $3678 


4 
x 
i 


9J bedicoebs snsbail 3? absA Ban inzesee (rae 


a @ 


- 
2 he 
o aa 


20 od 3 maol nois7z& asad toasass Sinspso | 


snos afd ni Bayet es¥ Ysivitos Bawwasem Je 


judixz2aib Gevyesdo Sat. é!.tii seri noitsoé 


(iW 


aviee bons eduioe Bedtcabse io WFhLide a 
pean 
¢ m 


teh sit .onidosasl siiiea fie Relzqieebs 2e.enneieas 
—_ 


wistsluriz2 sid o7 sonslanese $7 proige exaed. bedfieag 
me 
19732 io Jnemavom 362, ,G@e5 Li etuees ge 


a7 


‘(— —§#) TFos uozti0y 
Iq pue (@ —— @) [TOS uozTioyH dy weoy] uojJerIg 07 poaqiospe suUepUTT *SWAZYIOST SUPPUTT ‘HT°TI]T ean8tgq 
(qu/3n) uoTJerJUSsDUOD WNTAGTT I INby 
Oe edeak c 1 8°0 ed) 0 0 
a 
}- 
5 
Qu 
fo 
5 
()) 
© 
jan 
169) 
She 
o 
© 
AQ. 
ee 
e 
5 
H- 
rt 
= 
Gta. 
te) 
om 
ct 
10>) 
fo) 
h. 
oe 
OC ee 
— 
= 


Gc 


ye 


ova pe 


a0 


iso gory (8 —— @) 


r’e 


a ‘ 
¢ + i 
Ces RRS 0 OR, Ge A 6 ene 


2 . z = Be 


. 
rugume eqrokoeg\ mre soTSuc. sory (rE) = 


98 


0 
5 
e 
el 
fe 
4 
O. 
Y 
a 
i 
20 = 
0 . ae 4 le ae 1 
A/A (Activity / Total activity added) 
PiguretilalalS: sie activity distribution observed in May, 1982, after 
14 


C-lindane was applied to the Ap horizon (0 - 2.5 cm) in Dec., 1981. 
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adsorbed solute. 

69% of '*C-lindane activity applied at the upper 
Surface of the Ap horizon was accounted for by sample 
activities determined from the efficiency versus SCR curve. 
It was shown in chapter II that the activities calculated 
from the efficiency curve underestimated the tracer activity 
when counts were made in the soil sample-scintillation 
cocktail mixture. Volatilization would also account for some 
loss of activity (Kay and Elrick, 1967; Shearer ef al., 
1973). The relative contribution of each factor can be 
assessed by comparison with activities calculated for 
lindane added to the Bt horizon, which was covered in 15 to 
22 of moist soil, and would have suffered no volatilization 
losses. In the latter case 79% of the ayer Be hiew originally 
added was measured using the direct counting method. The 
difference between activities accounted for in the Ap and Bt 
horizons, 10%, suggest some volatilization losses from the 


surface layer. 


b. ‘'*C-lindane applied to Bt horizon 

The peak of mean activities of '*C-lindane added to the 
Bt horizon was also found at the original level of placement 
(figure F11.16)) Activity was aesecus five cm below and two 
cm above the level of placement. Amounts of '*C activity 
found below the level of placement may reflect problems in 
resolving the samples, which makes quantitative 
interpretation of these data difficult. The high apparent 
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14 


C-lindane was applied to the upper surface of the Bt horizon in 
Dec. 1981. Standard deviations shown. Depth coordinate given as 
distance above (negative values) or below (positive values) level 


of application. 
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from splitting of the layer of placement into two samples. 
Nevertheless, it can be concluded from these observations 


that lindane did not move appreciably in the Bt horizon. 


C. SIMULATION OF FIELD EXPERIMENTS 

A comparison between the simulation experiments and the 
field experiments is made in this section. Although the 
model was based on theoretical principles, it was built with 
the conditions of the field experiments specifically in 
mind. Therefore the set of observed experimental conditions 
does not provide a completely independent objective test of 
this model. 

The test being made is whether the model simulates the 
effects of this set of conditions on solute transport. The 
comparison may indicate what structural modifications to the 
model are required and what further experiments are 
necesSary, initially to provide information on how to adjust 
the model, and ultimately, to validate the model. 

To test whether the model will simulate the system it 
was intended to describe, it 1S necessary to use the same 
set of physical conditions as occurred in the field 
experiments. These conditions are discussed in the following 


Subsections. 


1. Physical Properties Controlling Simulation Experiments 


a. Bulk density 
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Dry bulk densities (table III.9) were calculated from 
wet bulk densities (determined from y-densitometry) and 
gravimetric moisture determinations, at 2.5 cm increments, 


as described in Chapter II. 


b. Porosity 

Porosities of Ap, Bt, and BC horizons were determined 
separately because of the structural and textural variation 
that occurs among them. Porosities, determined from 
pycnometer determinations of Dp and y-densitometer 
determinations of Db, were inconsistent with the moisture 


retention curve at many layers where a was calculated to be 
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at five kPa. Rawitz et aj. (1982) reported 
errors in y-densitometry resulting from use of the 
manufacturer's standards to construct a calibration curve. 
This is due to inaccuracy introduced by the assumptions that 
y absorption in soil can be defined by a universal average 
attenuation coefficient and that no backscattered radiation 
was measured by the detector. They reported serious 
ations of Db using the the same model two-probe 
density gauge as was used in this study. 

Overestimation of Db would result in underestimation of 
a when the equation a = 1 - Db/Dp is used. USing an average 
Dp value to calculate a within each horizon could also 
introduce error into calculations of a. Values of a which 
were calculated to be less than @ at five kPa are marked 


oo mind 4 ee =] 1 
with can’ (+) in*tablerEiIie9: 


' 
movi bosslusiss steme 
one (ezysmotiensh-_ 


—_ 
sf _ 


ainamersni mo 2.8 Ze , ane 


s1ujxes Gra lLesutousse oa3 Si on 


721 (S82) ie Fo asiwsll .694 ever ss be bem 


jeszabou ni tives? Bivow da ico sotsiamiseereee 


re ; 
aft dsiw tns3eleqnoonl sis# 8G to encktsan 


nieau .beeu 2a: gd\dad - | = » aoltseuos sole 
= Ad 


aan 
s19w snosizod DH Bos-,76 \GAq 


) 


rr 
ww 


1imteseb .eeliisemrd ands enon 


@- 


i9h-y Boe ga 36 enotsenimadds 


; =e - 
ois 26” » Stew etsysk-_asm 36 BVIUG Hoes 
§, 


ro 


4 
; y 7 f rt ry > ; mo 35 r ¥ 8: me 
A Aw " SJ no 4 ; {- i 4 J ic a 4a ° , akon n 


‘tien 6 rourtancs of eBbsabnadga ean 
26 
ea oft vd Bscyubotsnt yosesussent oF See. 
1 ore , 
s ¥d Seolte sd nso liow ni neisqaaeue 

>= 
ccestoed on t6f¢ bos snetolizess #ehseoum 


FP 


oliee betszoqges vernT .t06s0687eb sda vd beiavesen 


hom emee ait efig oniau dad do enolzemts 19 
-vbute efds nit Geaev @ew 26 Susp qi eT 


a) 


? 
> a 


io i 


sfyos sositod dose midntiw « stefusiss 6o2 Stee 


: ae 
2eulaVv .a to snetzéelusias o3at 16798 S203 
S 


$16 st svi? ge 8 neds anal 9d. oF bedalyales 5? 


e. lil side a be) 


Table II1I.9. Physical properties of Breton Loam. 
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To account for this discrepancy in the simulation of 


field experiments WGRAV was considered to be no less than 


2072 cm/layer (2028 -% bulk volume) 


in the Bt horizon and no 


less than .052 cm/layer (.020 X bulk volume) 


horizon. These values were estimated from the moisture 


in the ‘BC 


retention curve (figure III.17) as the difference between 86 


at 0 kPa and @ at 33 kPa tensions. 


c. Soil moisture retention curve 


Stagnant and mobile water pools were divided on the 


basis of soil moisture tension curves as discussed in 


Sections ll.A. Morstusesretention curves, (figure. bll.17) 
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were determined for Ap, Bt, and BC horizons. Celie te 
values were converted to volumetric by multiplying by the 
bulk densities given in table III.9. 

The stagnant water pool, WR, and the mobile pool 
Capacity, WMCAP, used in the simulation model, are heights 
of water (mm) corresponding to 6 values listed in table 
III.10. WR was determined as the product of 8 at either 200 
kPa or 500 kPa, and the layer thickness. WMCAP was 
determined as @ at 33 kPa moisture tension times layer 


thickness, less WR. 
Table III.10. 
Volumetric moisture contents used to determine WR and WMCAP. 


8 fat soil ftension of: 
HORI ZON DEPTH (cm) 33 kPa 200 kPa 500 kPa 


Ap Bare Dis) UR 366 0.290 Zee, 
Bt Oman c 0.418 Oso al 0. 32 
BC 60 0.384 Ooo Oi a0 


Gntiltration 

Initial conditions for the simulations of field 
experiments were based on 6 profiles measured on Mar. 9, 
Mar. 20, and’ Mar.25. Net infiltration occurred’ from. Mar 720 
to Apr. 26, while net evaporation occcurred from Apr: 26 


until the end of the experiment (table III.11). 
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Table III.11. Infiltration measured at experimental plots. 


Total change in Average daily Average daily 


Measurement water content infiltration evaporation 
period (mm) (mm/day ) (mm/day ) 
MareSea 20 Om00 = S 
Mem-20h- 25 (+)2505 0.41 = 
Mace2S5eApr 08 (490127 OFS = 
Apr.17 - 20 (+)7.50 205 - 
Apre20t -t238 (=): FG $.2 = 
Rie 2 3 3-826 COR S07 = 
Apr.26 - 29 (-)0.72 : 0.24 
Apr.29-May 3 Cael ere = 351 

May S2eta5 C=) 8 Gin? co #38 


The sensitivity of the model to the resolution of 
infiltration data was discussed in section III.A. During the 
Penlodsormmost intense infiltration, from Apr j917 to Apr. 
26, ® profiles were measured every three days. The relative 
daily rates of infiltration were determined by plotting the 
daily rates averaged over the period of each measurement, 
drawing a smooth curve through the points and interpolating 
as shown in figure 11.5. 

Runoff occurred under the snowpack at the 
y-densitometer access tubes. No runoff occurred at the 
points of application of the tracers because these sites 
were contained by heating duct extending above the soil 
surface. Infiltration where the tracers were applied was 
therefore greater than where infiltration was measured. To 
account for the difference the daily infiltration rates as 


determined by y-densitometry were increased by the factor of 
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the volume of water determined by the Mar. 25 snow survey 
divided by the total volume of infiltrating water measured 
at the access tubes. 

Snow drifts along the east side of the plot area 
covered some sites until Apr. 29. As much as two times as 
much water was held in these drifts as occurred on the west 
side of the plots. The soil under this snow also thawed more 
Slowly. The y-densitometer access tubes on the east side of 
the plots were kinked during installation and so the data 
necessary to describe leaching here were lacking. Therefore 
only results obtained for the experiments on the west side 
of the plots, where a complete set of simulation input 
variables was available, were compared with the results of 


Simulations. 


e. Frozen layers 

The depth of thawing was determined by physical 
probing. However this method did not indicate the maximum 
depth of thaw, which fluctuated dvurndliy. ~Thetdepth to 
frost measured by probing did, however, coincide 
approximately with the bulges observed in the differences 
between successive 6 profiles. Gillies (Gray et al., 1970) 
observed that during infiltration of Prairie snowpack two 
distinct layers developed: an saturated thawed layer 
overlying an unsaturated frozen layer. The bulges between 
Successive 6 profiles was used to approximate the maximum 


depth of daily thaw in simulations of field experiments. 
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2.Comparison of Simulations and Field Observations 

Simulations were performed which paralleled the field 
experiments described in Section III.B using the measured 
and derived parameters and variables described above. Data 
files for two of the simulation runs (1.b, *HOH applied to 
Phe yBurRhorizon: andmzvar “£C=lindane applied to the Ap 
horizon) are given in Appendix E, along with the final 
output generated by running these files in the simulation 
program. The FORTRAN program which implements the simulation 
model is also printed in Appendix F, with instructions for 
use in Appendix D. 

Observed data used for comparison are expressed as 
A/A,, where A is the activity in each layer and A, is the 
total activity detected in the soil column. The simulation 
model is concerned only with vertical transport of solute 
and relative activities are a truer measure of vertical 
displacement than absolute activities which reflect both 
vertical and lateral displacements. 

This convention has the additional advantage of 
correcting for the error introduced by counting activities 
directly in the soil sample (as described in chapter II). 
Accountability experiments showed that measured activities 
were proportional to known activities for standard soil 
samples (figures II1.2 and II.3). Therefore: 


A/A, = A*/Ao, 


where A measured activity in sample, 


A, total measured activity in soil column, 
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Ax 


actual activity in sample, 


Ao 


total actual activity added to soil column. 

The data are reported as mean A/A, values for each 
layer in the graphs which follow and in the corresponding 
tables in gAppendix, CC, except for figure II1.18, for which 


only one set of data was obtained. 


Test of Fitness 

The agreement between observed and simulated data is 
assessed by regression analysis. The fitness test is whether 
the regressions of observed on predicted data have 
intercepts which are close to zero and slopes close to unity 
(Hillel, 1977). The regression correlation coefficients (r?) 
give an indication of the variation in the data from the 
regression line. Regression analyses were performed on a 
point-by-point basis, using all samples with measureable 
activity (100 total counts over background) and all 


predicted activity greater than 0.1% of A,. 


Experiment 1: Leaching of a non-reactive solute (°HOH) 


a. °HOH applied to Ap horizon 

Less than 0.3% of the tritiated water added to the soil 
surface in November and December, 1981, was accounted for in 
the samples taken. Simulation of the experimental system 
with spring evaporation only did not predict this loss. 
Though drying of the soil surface was observed in early 
December, no evaporation data were collected at this time. 


Early winter evaporation was simulated and a better 
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Figure III.18. Observed ( ) and predicted (— —~) spring 
tritium activity distribution for 3H0H applied to Ap horizon (0 - 2.5 


cm) in the fall, for the one site where activities were detected. 
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prediction was obtained, though evaporation data for this 
Simulation was speculative. Simulating daily evaporation of 
0.1 cm over ten days in December predicted that 76% of the 
tritium activity would be lost before infiltration began in 
the spring. 

Two activity peaks occur in the observed tritium 
activity distribution (stepwise plot, figure III.18). The 
lower peak iS approximately predicted by the simulation 
(smooth curve). The depth of penetration was also reasonably 
predicted by the simulation. The possibility that the upper 
peak was caused by microbial assimilation was discussed in 
the previous section. If this peak activity is disregarded 
then a much better fit between observed and simulated solute 
distribution is seen. Regression of observed distribution on 
predicted is given below. 

Including upper peak: 
Obsoes 0471 +8. 174P Sr 2ue?d .04% 
Excluding upper peak: 


Obs 225ni0 15S +8 2159% Pom rbe=oX 62: 


b. °*HOH applied to Bt horizon 

Three different conditions were simulated for 
comparison with observed distribution of tritium activity 
for *HOH added to the Bt horizon: A) using 200 kPa soil 
moisture tension to assign WR and WMCAP; B) using 500 kPa 
tension to assign WR and WMCAP; C) disregarding the effects 
of a frozen layer. One set of the data used in simulation A, 


along with the output generated, iS given in Appendix D. For 
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B, WR and WMCAP were changed from the values used in A 
according (fo Onvaluesrqivenvin tables J11.10.\ ror :C internal 
variables were identical to those used in A, but FLNR was 
set at zero at all time steps. All three simulations were 
subject to the same infiltration inputs. 

The results of the simulations were compared with the 
results of the field experiments performed on the west half 
of the plots (figure III.19). Observed activity distribution 
at one site indicated that net upward transport of solute 
had occurred. This data did not fit with the observations at 
all other sites, where net downward transport occurred. It 
was concluded that some mechanism other than those being 
modelled was operating at this site and this data was not 
used in the comparison with the simulation resuliees 

Comparison of simulations A and B with field 
observations (figure III.19) showed a qualitative agreement 
in shape and position of the solute activity profiles. All 
three solute profiles possess a high degree of symmetry 
about a point seven to eight cm below the level of tracer 
application. Simulation A gives a better prediction of the 
position and value of the mean activity peak. Both A and B 
accurately predict the depth of solute penetration. 
Regression of the observed tritium activities on those 
predicted by the model are given below: 

A (WR:WMCAP at 200 kPa tension): 
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B (WR:WMCAP at 500 kPa tension): 
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Observed 


Predicted: 


—e ee eo BOA MT ba Ona 
(WR:WMCAP set at 2 bars tension) 


Simulation B 
(WR:WMCAP set at 5 bars tension) 


ee eo LIN ba COT EG 
(no frozen layer). 


Depth (cm) 


a) 503 220 BES, 220 ap ae 


A/A, (Activity / Total activity present) 


Figure III.19. Mean observed and predicted spring tritium activity 
distribution for ROH Aperidd to upper surface of Bt horizon in fall. 
Depth coordinate given as distance above (negative values) or below 


(positive values) level of application. 
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Obsi= D022 20. OOxP ort ase, 730% 
Both regressions have intercepts close to zero and the 
correlation coefficients are nearly the same. The slope of 
the regression line was closer to 1 for A, indicating a 
Slightly better fit with the observed results. 

‘Comparison of the use of 200 kPa and 500 kPa soil 
moisture tensions to assign the sizes of the mobile and 
Stagnant pools was inconclusive, mainly because of the 
influence of the frozen layer discussed below. The 200 kPa 
value produced a marginally better fit with the observed 
data and is used for simulations of transport of the 
adsorbed solute. 

Neither simulation A nor B predicted the presence of 
solute more than five cm above the zone of application, but 
5.5% of the activity detected in the soil was found above 
this depth. The possible causes of the upward displacement 
are diffusion, evaporation and thermal gradients. 

Mixing in the model is achieved by layer thickness and 
occurs only when there is movement of the soil solution. 
Total spreading effects were shown in section III.A to be 
relatively insensitive to layer thicknesses between 1.3 and 
ten cm so that increasing Az up to ten cm would not cause a 
large increase in mixing (see figure III.7). Changing Az to 
effect greater diffusional mixing would also spread out the 
concentration profile at the advancing front where a 
reasonably good fit is already achieved. To model diffusion 


more accurately, for both moving and stationary solution, 
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the model would have to incorporate a specific rate law 
describing diffusion (Fick's law), along with some means to 
increase diffusivity (D) when moisture content is high, such 
as the approximation formula presented by Shearer et al. 
(1973) given in the literature review. 

The moisture potential gradient produced by evaporation 
was not measured but large evaporative losses of water were 
observed in 6 profiles. Temperature gradients which occurred 
in early winter would have been small (soil water 
temperatures generally vary between 0°C and 5°C at this 
time). During spring thaw temperature gradients are opposite 
in sign and would result in downward moisture flux. 

No changes should be made to the model to account for 
the effects of diffusion or evaporation until experiments to 
isolate their influences have been conducted. 

The symmetry of these solute profiles emphasizes the 
influence of the frozen layer on solute transport. The 
Simulated effects of structure, namely an upwardly skewed 
solute profile (simulation C), are not evident in the 
Simulated nor in the observed results for the partially 
frozen structured soil system. The frozen layer acts as a 
control on infiltration rates causing more complete 
equilibration of solution between mobile and stagnant pools, 
and producing the more symmetrical solute profiles seen in A 
and B. 

Comparison of A and C, which possess the same internal 


properties except for the frozen layer, show that it is 
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possible for solute transport to be greater for leaching of 
soll-borne solutes during thawing of a partially frozen 
Structured soil than for leaching of an unfrozen structured 
soil. This effect would occur because the infiltrating water 
is held in contact with the intra-aggregate solution by the 
impeding frozen layer thus picking up a greater solute load 
from the soil. Majka and Lavy (1977) described a similar 
mechanism to explain observed increased solute transport in 
packed columns of fine textured soil over coarse soil. They 
Stated@thateinethisecasemtexeuresexerts a0limitingsinfituence 
on flow rates thereby causing increased equilibration 
between soil-borne adsorbed solute and solution phase. 

The above comparison (A versus C) suggests another 
possible explanation for the observed solute distribution. 
During Slow infiltration of meltwater clay soils are apt to 
swellp( the tBevhorizontoftthisssoil r6%33%telay)varn “such 
cases the intra-aggregate macropore volume may be sealed off 
by the expanding soil matrix. Hydraulic conductivity and 
flow rate are reduced; equilibration 1s more complete and 
the resulting activity profile is more symmetrical as shown 
iis ne Syme Lnymormine sObServed SOlULesCISEerImuLTOn could 
result from leaching through structure or non-structured 
partially frozen soil. Because the effects of the frozen 
layer mask the effects of soil structure it is necessary to 
study the effects of structure in isolation before definite 
conclusions can be made about the the influence of structure 


alone on solute transport in systems with well developed Bt 
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horizons. 

Using a solution to the Darcy equation, rather than a 
water balance equation, to determine flux could possibly 
have determined the mechanism which caused the symmetrical 
solute concentration profile observed. This would have 
necessitated calculation of both saturated and unsaturated 
hydraulic conductivities (K) that vary with swelling of the 
Soil. Because of the heterogeneity of the soil this would be 
dvtiiculestoyaccomplishn in the field) Laboratory 
measurements could be used to measure K for each horizon, 


but are much less reliable than field measurements. 
Experiment 2: Leaching of adsorbed solute ('‘C-lindane) 


a. '*C-lindane added to Ap horizon 

Less han 2pOte aca lindane ;ACGEG ptOutnew Op wo. cm Or 
the Ap horizon was detected as '‘C activity beyond this 
layer (figure III.20). Linear regression analysis of 
observed data on predicted data generated the following 
regression equation: | 


ODSh=) 204 see ato eR r= OO92:, 


Eliminating the effects of soil structure from the 
Simulation (that is setting WR=0) predicted the observed 
lindane distribution almost equally as well, however 

(Obs =0-.052 HlwlSosPoe tg 9. 98/).¢ these results «point out 
the importance of adsorption on solute transport. A 


moderately high adsorption coefficient (12.6 ml/g for Breton 
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Figure III.20. Mean observed (——) and predicted (— —) spring 
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Loam Ap horizon) overrides the effects of soil structure in 
the simulation experiment (see also figure III.9). The 
dominant influence of adsorption is also reflected in the 


observations of the field experiment. 


b. '*C-lindane added to Bt horizon 

The data obtained from distribution of ‘'‘*C activity 
resulting from transport of labelled lindane added to the Bt 
horizon were less clear cut, mainly because of the 
difficulties in resolving the field samples. At two of the 
Sites used for comparison with simulations the original 
layer of placement of lindane may have been split between 2 
samples resulting in a mean distribution which may overstate 
the actual solute movement (stepwise solid line plot, figure 
III.21). The Simulation (smooth curve) predicted that more 
than 80% of applied lindane would be found in the layer 
where it was placed. The uncertainty in sample resolution 
makes interpretation of this experiment difficult, however, 
it can be seen that adsorption continues to dominated the 
distribution of solute in the Bt horizon under the leaching 
conditions studied. 

Beyond restating the importance of adsorption in solute 
epCaaeet little need be said about the experiments with 
'4C-labelled lindane. Simpler models predicted the observed 
distribution of solute as well as the structured soil-frozen 
layer model presented here. Experiments performed under more 
intense leaching condition are necessary to assess the 


effects of soil structure on transport of adsorbed solutes. 
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Figure III.21. Mean observed ( 
oc activity distribution for eee dane applied to upper surface 
of Bt horizon in the fall. Depth coordinate given as distance above 


(negative values) or below (positive values) level of application. 
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Validation testing 

The experiments discussed above indicate that the 
mechanisms modelled in the program require further 
elucidation before complete structural validation of the 
model is possible. Experiments designed to study each 
mechanism in isolation are required. Simulation experiments 
can prove useful in pointing out the experimental conditions 
which are likely to provide information necessary to test 
the model. That iS, variation of input which causes definite 
differences in simulations can be used in physical 
experiments. The observed results of the physical 
experiments will then permit conclusions about the 
structurale validatyeotetheesimulatiaonqamodely orl partsnofait. 

Once the important mechanisms have been elucidated and 
included in the model research should be performed in the 


field to test the predictive validity of the model. 


Tat 


ed+ jad? ss65%bal 
ted3qi e1iupe’ BIPOT oi + 

ens 3o moivebilev isiutserge sdelq 

doses ybuia oF benples® aaneatiisax 
a¢nemitsqxes nOlijelumi2 -bsyinpe? sre 6! 
tin [sjnemisegxe sit sue patsarog ni Ig 
{is22ecsn folseqrotal abivoig of 

siinitaB esses doldw svaqnt Ja nelsaizaw ai aes 
vig of Beev sd.iso soolgsliuela he 


eptevda ad’ te a+iuess i ap sit 22 


' 


utmios jinizas nent Ltiv a aan? 

e 3 fisq 20 .lebom nolitslumis sds 36 7 tbe daw Bs 
324 svat sm2iocsiset tist200mi dai 
fe dotee@et [sbom agg wile 


shom scz to viibilev svidoibew stg a on 


: 
7 a 
aN 738 


IV. SUMMARY AND CONCLUSIONS 


A. SUMMARY 
Model Description 

A finite difference simulation model based on the 
chromatography/continuity equation has been developed to 
describe solute transport through structured soil, anda 
"FORTRAN' computer program written to implement the model. 
The program can be used to Simulate transport of either 
adsorbed or non-reactive solutes through partially frozen or 
completely thawed soil. 

The model characterizes a heterogeneous soil solution 
and pore geometry by defining two interactive solution 
pools: mobile and stagnant, based on the concepts presented 
by Addiscott (1977). At each infiltration event only the 
mobile solution is displaced. After each displacement a new 
equilibrium is established between the pools within each 
layer. 

Transfer of solution between layers is defined by a 
water balance equation adapted from Burns (1974) to allow 
for variable water content. The water balance equation is 
extended to describe the influence of a frozen layer which 
blocks transmission of water and causes water storage above 
the frozen layer in excess of field capacity. A first 
attempt to model the effects of evaporation, based on the 


method used by Burns (1974), has been incorporated into the 
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model. This method was extended to describe the mechanism of 
volatilization of solute from the surface layer of the soil 
system. 

Transport of adsorbed solutes is treated by partition 
of solute between solution and solid phase on the basis of 
an effective partition coefficient determined from a linear 
adsorption coefficient, bulk density, and volumetric 
moisture content. 

The computer program which implements the model is 
reproduced in Appendix F, along with a users' guide and 
Sample input and output data (Appendices D and 5, 
respectively). 

Simulation Experiments 

Simulation experiments were performed to assess the 
sensitivity of the model to the relative sizes of mobile and 
Stagnant solution pools; infiltration rates and pattern; 
layer thickness; frozen layers; and adsorption of solute. It 
was found that the model was sensitive to adsorption 
coefficients, infiltration rates, and presence of a frozen 
layer; but not sensitive to layer thickness variation 
between 1.3 cm and ten cm. 

Simulation of Field Experiments 

Transport of tritiated water (non-reactive solute) and 
'4C-labelled lindane (adsorbed solute) through partially 
frozen soil during infiltration of snow melt was studied in 
field experiments performed on a Gray Luvisol. Simulations 


of the field experiments were run using the data obtained 
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from analysis of the soil physical and chemical properties. 

Agreement between predicted and observed results were 
assessed by regression analysis. The simulation predicted 
the peak concentrations and downward extent of solute 
movement reasonably well for both adsorbed and non-reactive 
solutes. It did not predict the upward extent from the zone 
of application of the non-reactive solute. This lack of 
agreement indicates that the mechanisms of upward movement 
caused by evaporation and/or diffusion should be elucidated 
and incorporated into the model. 
Future Experiments 

Conclusions drawn from comparison of simulations and 
field observations regarding specific mechanisms controlling 
solute transport are tentative because of interactions 
between the factors involved (soil structure, frozen layer, 
infiltration rates, adsorption-desorption). Some of these 
factors produce opposite effects on solute transport; soil 
aggregation causes skewed solute concentration profiles 
while frozen soil causes symmetrical profiles. Other factors 
reinforce each other; adsorption and aggregation both cause 
holdback of solute. Experiments which study each factor in 
isolation are necessary to validate the treatment of each 
mechanism by the model and to test the validity of 
assumptions made in constructing the model. Simulation 
experiments can be used in this context to indicate the 
experimental conditions which would best elucidate the 
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predictive leads (Hillel, 1977). 

The assignment of values to define the sizes of the 
mobile and stagnant solution pools is fundamental to this 
and other models (Addiscott, 1977) of solute transport 
through structured soil. Optimization experiments using 
undisturbed soil columns are necessary to determine the 
appropriate soil moisture tension on which to base this 
division. The model makes the assumption that equilibrium 
between solution pools occurs after each infiltration event. 
Experiments designed to test the validity of this assumption 
under various moisture flow rates are also necessary. 

After the factors and mechanisms influencing solute 
transport have been adequately defined and mathematical 
expressions of their effects included in the model the 
experimental venue should be moved back into the field to 


test the predictive validity of the model. 


B. CONCLUSIONS 
Simulation Experiments 

The simulations performed support the following 
‘conclusions: 
i) Solute held within soil aggregates are less subject: to 
leaching by infiltrating water than are soil-borne solutes 
in non-structured soil. Solute distribution in structured 
soils is characterized by an assymetrical concentration 
profiles (upward skewness) while solute concentration 


profiles in non-structured soils are symmetrical about the 
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average displacement of the solute. 

1i)Adsorption plays a dominant role in controlling solute 
transport through soil, even when adsorptive forces are 
relatively weak. 

iii) The rate of infiltration influences the leaching 
pattern which develops in structured soils to a greater 
degree than in non-structured soils. Intensive infiltration 
transports less soil-borne solute but carries a small amount 
to greater depths than low intensity infiltration. This is 
reflected in skewed solute profiles. 

iv) A frozen layer influences solute transport by slowing 
infiltration and causing more complete equilibration between 
solution held within aggregates and solution moving between 
aggregates resulting in symmetrical solute profiles. This 
can cause transport of greater solute mass under some 
conditions of leaching in highly stuctured soil. The slow 
infiltration of meltwater into partially frozen ground may 
represent the most intensive leaching conditions encountered 
in some structured soils. 


Field Experiments 


i) Non-reactive solute: 

The simulation of transport of tritiated water through 
partially frozen, structured soil agreed with the observed 
distribution of solute from the zone of solute application 
downward, in peak concentration position; slope; and maximum 


depth of penetration. Simulations did not adequately 
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describe the solute distribution above the zone of 
application, indicating the mechanisms of upward solute 


transport of solute are not adequately defined in the model. 


ii) Adsorbed Solute: 

Adsorption coefficients of 12.6 and 5.1 ml/g were high 
enough to exert a dominant controlling influence on 
transport of ‘*C-lindane through the Ap and Bt horizons, 
respectively, for the leaching conditions observed during 
infiltration of meltwater into a Gray Luvisol. Predicted and 
observed results both supported this conclusion with a high 


degree of correlation. 
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APPENDIX A 


List of Symbols Used in Chromatographic Equations 


N 
Hl 


vyerticatieds stancel coordunateitcm) . 


> 
" 


e@ross=sectional areaelotdson lucotumny (em? )& 
Me=ivolume,~AAz cof TayervofathieknesstAz) (cm*). 


Vp = pore volume, a-V, of layer of thickness Az (cm?). 


pest timen(d). 

eS concentrationuwa7m)) . 

pee water flux densivy (Darcian, velocityIé(ml/cm7-d) . 
y = water velocity (cm/d). 


of] solute mass. flux@density (g/cm’-d). 


eolnee ditfusivanry in soil ((cem?/d). 
Dee="SOLUte ditiiusiyity in water (cm?7d). 
EP dispersion coetficient (cm’/d). 


Dspr = spreading coefficient, combining D and E (cm?/d). 


T = tortuosity factor (dimensionless). 

Ge= Porosity (cmG/cm-).. 

6 = volumetric water content (ml/cm?). 

6s = saturated volumetric water content (ml/cm°). 
Db =Adry bulki densitytqvemen . 


a = mass adsorbate/mass soil (ug/g). 
Kads = adsorption coefficient, Freundlich isotherm (ml/g). 
B = effective partition coefficient (dimensionless). 
=e Db ky 6. 
The symbolism used in the simulation program 1s given 


in Appendix D. 
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APPENDIX B 


Summary of Field Data 


Table B.1. Mean radioisotope activity distributions observed 
in May, 1982, after labelled solutes applied to the Ap 


hemaizonre 0se1t23) cm) in) Nov. and Dec.:, 1981. 


Depth (cm) A/Ao(?HOH) ' A/A, (' *C-lindane) 
OTS 23 O75 O78 ty 206 
Shee | Cat .0049 OMA 10 2.6 
St ee 0 Oye ee 
ier 0 «0029 
12a 125 .0029 

2 3et 5 0031 
im. -— 11s - 0.025 

1s = 20 0 

AO 22 25 oi LOZIS/, 

Lae ae 0 


A = activity measured in each soil sample layer. 


Ao = total activity added to each site 


6.25 uCi iow 4HOH. 


2 nuCi £Eormey2C—) indane:. 


‘Note - activities detected at one site only. 
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Table B.2. Mean radioisotope activity distributions observed 
in May, 1982, after labelled solutes applied to the upper 
surface of Bt horizon in Nov. and Dec., 1981. Depths are 
given as distances from the centre of the layer of 
application (negative values above layer of application, 


positive below). 


Depth (cm) A/Aj, (?HOH) A/Re * *C=bindane) 
mel ts ~UO2t 25.0041 

me p23 a0026 “EOU30043 

=r eae 10042) ta. 0067 

74 .0064 + .0094 

mas “0094; to). 0:18 
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20 a0) 257 Peet 227 
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25 s0004 + 3.0007 

215 a) O10 Aa 0.0.07 

A activity measured in each soil sample layer. 


total activity added to each site 
(6.25) anda 1t0muGaetor “HOH 2 1C1 (for “*C-lindane) 
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APPENDIX C 


Table C.1. Mean radioisotope activity distributions observed 
on west half of experimental plot area in May, 1982, after 
labelled solutes applied to Ap horizon (0 - 2% cm) in Nov. 


anadtDecen/ 1981. 


Depth (cm) A/At (*HOH) ' A/At GY *C= lindane) 
Uae cmer o2o2 OL eters ae 
2h OS oS .029 + .044 
Shae 3 0 One 0 
ice ae ie liz 
PO =e 2s lee 
ie eee hS 124 
(oa wales 099 
Piss 0 
20. 225 068 
Pee) at PRs 0 
A = activity measured in each soil sample layer. 


At 


total activity measured at each site 
= ZA. 


‘Note - activities detected at one site only. 
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Table C.2. Mean radioisotope activity distributions observed 
on west half of plot area in May, 1982, after labelled 
solutes applied to the upper surface of Bt horizon in Nov. 
and Dec., 1981. Depths are given as distances from the 
centre of the layer of application (negative values above 


layer of application, positive below). 


Depth (cm) A/At (?HOH) A/At(' *C-lindane) 
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APPENDIX D 


DIMENSIONING: 

The program is set up for systems of less than 100 
layers. If the system to be modelled contains 100 or greater 
layers the dimensions of the following arrays must be 
changed. 

WT(m), WR(m), WMCAP(m), WGRCAP(m), WM(m), WGRAV(m), 

XSWMCP(m), XSWGCP(m), ST(m), SR(m), SM(m), SGRAV(m), 

WE(m), SE(m), MS(m), Q(m), DB(m), KADS(m), DEPTH(m) 


where m = number of layers in system. 


PROBLEM IDENTIFIERS: 
TITLE1 = first title line, up to 80 alphanumeric characters 


TITLE2 = second title line, up to 80 alphanumeric characters 


SYSTEM PARAMETERS: 
M = number of layers in system being modelled. 


DELTAZ = layer thickness, in mm. 


TIME STEP PARAMETERS: 
TFIN = time period of simulation, days; 

also signals when the final output data 1s printed. 
T1, T2, T3 = times when intermediate results can be printed; 


if not required set all equal to TFIN. 


LAYER PARAMETERS: 


DB(m) = dry bulk density of each layer (g/cm*); used 
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only for adsorbed solutes, otherwise can be set = 0. 
KADS(m) = adsorption coefficient, linear Freundlich isotherm 
(ml/g); must be specified (for non-adsorbed solute 


set <=) 0). 


INTERNAL VARIABLES: 
Initial values for the following variables must be 
provided by the user. 
MS(m) = total mass of solute in each layer per unit area 
(e,g., kg/ha, ug/em*): area can be chosen 
arbitrarily, but must be consistent with units 


used for S (see DRIVING VARIABLES). 


WT(m) = total height equivalent (mm) of water in each layer; 
determined by user as product of volumetric moisture 
content and layer thickness in mm. 
= DELTAZ*® 
WR(m) = height equivalent (mm) of water in stagnant water 


pool; determined from moisture retention curve (8 at 
200 kPa moisture tension) and layer thickness in mm 
= DELTAZ*6(200 kPa) 
WMCAP(m) = mobile water pool capacity, in height equivalent 
(mm) of water; determined from moisture retention 
curve and layer thickness in mm. 


DELTAZ*[86(33 kPa) - 86(200 kPa) ] 


WGRCAP(m) 


gravity water capacity in height equivalent (mm) 
of water; determined from porosity, moisture 
retention curve and layer thickness in mm. 
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layer routine only; can be set at 0 at all 
layers for non-frozen soil. 

The simulation program generates the following initial 
internal variables from those provided above and updates 
them throughout the program: 

WM(m) = water contained in mobile pool, height equivalent 
(mm). 
WGRAV(m) = water contained in gravity pool, height 
equivalent (mm). 
XSWMCP(m) = unfilled capacity in mobile water pool (mm). 
XSWGCP(m) = unfilled capacity in gravity water pool (mm). 
B = effective partition coefficient (dimensionless); 
= DELTAZ*DB(m)*KADS(m)/WT(m). 


Q(m) = mass of adsorbed solute (ug) in each layer; 


MS (m)*B/(1+B). 
ST(m) = mass of dissolved solute (ug) in each layer; 
= MS(m)/(1+B). 


SR(m) = mass of dissolved solute (ug) in stagnant pool; 


ST(m)*WR(m)/WT(m). 
SM(m) = mass of dissolved solute (ug) in mobile pool; 


ST(m)*WM(m) /WT(m). 


SGRAV(m) = mass of dissolved solute (ug) in gravity pool; 


ST(m)*WGRAV(m) /WT(m). 


DRIVING VARIABLES 
The following input variables must be provided at each 
time step. There must be TFIN sets of input data. 


FLNR = frozen layer number (dimensionless); 
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the hippersMiimms roof ftrost in the soil; 
for unfrozen soil must be set = 0. 
Wil=<darly ‘preciipation or infiltration (mm); 
(see note’, below). 
E = daily evaporation (mm) (see note', below). 
S = daily mass of water borne solute applied 
per unit area of soil; must be compatible 
wath units used for MS(m) “(edge , kg/ha, ug/cm7). 
' - Note: Net infiltration is determined from W and E as 
WA = W - E. WA cannot be 0 because some concentrations are 
determined by dividing by WA. Therefore W and E must not be 
equal. 

The following driving variables are generated by the 
Simulation program to affect transfer of water and solute 
mass between layers: 

WA, WL = water added from layer above and lost to layer 
below, respectively, by downward transmission of 
water infiltrating at surface layer. 


Sao 


solute mass carried by WA and WL, respectively. 

WE(m) = water lost to layer above due to capillary flow 
induced by evaporation; (WE(1)=water lost from 
surface layer to atmosphere by evaporation). 

SE(m) = solute mass carried by WE(m); (SE(1) = solute lost 


from surface layer to atmosphere by volatilization). 
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UNePUT CARD) INSTRUCTIONS. 
The data enterred into to the program must follow the 


following tformat: 


CARD VARIABLE FORMAT COLUMNS 
| TITLE! 80A1 Us ae ae ke, 
2 TITLE2 80A1 Pees) 
3 M LS bined 

TEIN is Gaaw hO 
sey LS Dit FES 
ae ES Shy ae a8 
ee is Fed ie PR 
DELTAZ BS 2 Zones 0 
4? DB P0445) 1 - 10 
KADS F10.4 Of >- M20 
MS F10.4 Ze es 0 
WT F10.4 ay 40 
WR F10.4 fier 50 
WMCAP F10.4 S67 e60 
WGRCAP F10.4 ati Pee eae Ae 
oe FLNR ins Loe he 
Ga W ELOw2 ee el 
E EERO £2 ier? 20 
Ss) FspOr a2 Zi oO 
Notes: *7 - Card 4 is repeated M times, where M is the number 


of layers in the system being modelled, in order 
that variable values are supplied for all layers. 


> - Cards 5 and 6 are repeated TFIN times, 
where TFIN is the number of time steps simulated. 
Cards 5 and 6 must be enterred alternately, so 
that daily FLNR, and W, E, and S are read 
together. 
Examples of complete input files illustrating the 
correct formatting can be found in Appendix E. Input data 


are enterred on 1/0 device unit ‘5; output data can be 
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OUTPUT FORMAT 

The results of the simulation program are printed in 
the following format: 
1. Model identifier (2 lines). 
2. Problem identifiers, reprinted from input file (2 lines). 
3. Time step and layer parameters, reprinted from input file 
(1 line). 
4, Data echo: the data for each layer is reprinted as it 
appears in the data file, under the headings, DB, KADS, MS, 
WT, WR, WMCAP, WGRCAP, as defined above. 
5. Initial conditions, calculated from input data supplied, 
under the headings, I, DEPTH, XSGMCP, XSWMCP, B, Q, ST, and 
SM, as defined above. 
6. Infiltration/evaporation Se mea (daily), under the 
headings J (time), FLNR (frozen layer number), FLNRT 
(previous frozen layer number), W-INPUT 
(precipitation/infiltration), EVAPN (evaporation), S-INPUT 
(solute carried by infiltrating water), W-POND (height of 
ponded water), WA(NETW) (net infiltration/evaporation), and 
WIOTAL (total imtiitrataon-tordatre), in tabular ttorm. 
7. If there is any volatilization loss of solute the above 
table is interrrupted and the following information is 
supplied: NET E (always 0); WT(1), WG(1), WM(1), WR(1) 
(total, gravity, mobile, and stagnant water, respectively, 
contained in layer 1); WE(1), SE(1) (evaporative water 
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Aye 
8. Sink losses: water and solute which drain from the bottom 
layer of the system will also interrupt the data tabulated 
in 7 (above) under the headings: SINK LOSSES AFTER DAY "J"; 
VOLUME, AMT. SOLUTE, and CONC. SOLUTE (drainage water, 
height equivalent; mass of solute drained, per unit area; 
and solute concentration in drainage, respectively). 
8. The solute concentration profile is printed in tabular 
form under the headings: LAYR (layer number), GRAVHOH, 
MOBHOH, and RETHOH (gravity, mobile, and stagnant water, 
respectively, contained in the layer): GRAVSOL, MOBSOL, and 
RETSOL (solute mass, per layer, contained in gravity pool, 
mobile pool, and stagnant pool, respectively); DISSOL, 
ADSOL, TTLSOL (total dissolved solute, adsorbed solute, and 
total solute mass, respectively, contained in the layer). 
Examples of the output data file are found in Appendix 


E, following their respective input data files. 
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APPENDIX E 


Input Data Used by, and Output Data Generated by Simulation 
Model Program. 

The following pages contain I/O data files 
used/generated by the FORTRAN program which implements the 
Simulation model. (The program is printed in Appendix F.) 

Simulation E1, pages 146 to 149, is for transport of a 
non-reactive solute (tritiated water) applied to the Bt 
horizon and is a simulation of field experiment 1b, 
dGescrmbed In@chapter oul). Samulation:.B2, pagesj150 to. 753, 
is for transport of an adsorbed solute ('*C-lindane) applied 
to thesgAp hoprzon (005 2.5 em)2 and sta simulation of field 
experiment 2a, chapter-III. 

The output file generated by each simulation follows 
its respective input data file. Details on formatting of 


input and output files are given in Appendix D. 
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Nel 


009 
4,20 
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149 
OUTPUT DATA’ FILE El; 


PFSS2: TWO-POOL MODEL OF SOLUTE TRANSPORT THROUGH 
PARPLALLYS FROZEN<STRUGTUREDESOIL. 


SIMUEATIONLOFORIELD EXPT. ib: SHOH APPLIEDITO Bt HORIZON. 
WREWMCAP \SETJAT 2° BARSSULAYERSH=140 X 25.4 imm. SITE.23. 

200. 5417.80 (400,050 25.4 
(NOTE:DATA ECHO, INITIAL CONDITIONS, INFILTRATION RECORDS, 
INTERMDIATE SOLUTE DISTRIBUTION PROFILES, ARE OMMITTED) 


SOLUTE PROFILEY DAY 54 
LAYR GRAVHOH MOBHOH RETHOH GRAVSOL MOBSOL RETSOL DISSOL ADSOL TTLSOL 


1 O20 0.0 728 0.0 0.0 B20 0.0 0.0 0.0 
2 0.0 0.0 138 0.0 0.0 0.00 0.00 0.0 0.00 
3 0.0 0.0 FSS 0.0 0.0 OrO1 0.01 0.0 Oeren! 
4 0.0 0.0 738 OO 0.0 0.08 0.08 0.0 0.08 
5 CFO 0.0 7738 0.0 0.0 0.40 0.40 0.0 0.40 
6 0.0 a0 7.938 020 0.0 £28 LeZs 0.0 1.26 
a 0.0 0.0 9.42 0.0 0.0 3.93 3693 0.0 3.93 
8 0.0 0.0 9.42 O20 0.0 8.38 8.38 0.0 8.38 
9 0.0 fb. 05 9.42 0.0 1,499 .+43.37 14.86 0.0 14.86 
10 0.0 Peel 9.42 0.0 2 0996. 26) leas 0.0 £8.35 
Vt 0.0 Pa2e 9.42 0.0 2. O039e93.84 “17788 0.0 17.88 
2 0.0 Ha 9.42 0.0 1.64 12.74 14.37 See 14.37 
13 0.0 Peli 9.42 0.0 Pelt O- 067 oy 0.0 ere 
14 C0 ew 9.42 0.0 0.65 an07 Dawe 0.0 Du he 
15 0.0 PoZd 9.42 O%O GeSo 2.58 2.91 0.0 2291 
16 Ox0 ez 9.42 0.0 Oka ie: 4D 1.30 0.0 Hisshe 
17 0.0 eZ. 9.42 0.0 0.06 0.45 0.59 0.0 Oeod 
18 0.0 aan! 9.42 0.0 0.02 0.167 “O..08 0.0 O21 
19 0.0 2 9.42 0.0 O01 0.05 0.06 0.0 0.06 
20 O20 Pet 9.42 0.0 0.00 opm en 0.02 020 02.02 
yaa O70 W221 9.42 0.0 0.00 0.00 0.00 0.0 0.00 
Za 0.0 ayaa 9.42 0.0 0.00 0.00 0.00 0.0 0.00 
23 O70 w.21 9.42 0.0 0.00 O00 0.00 0.0 0.00 
24 0.0 eZ 9.42 0.0 0.00 0.00 0.00 0.0 0.00 
25 O20 1.34 8.40 0.0 0.00 0.00 0.00 0.0 0.00 
26 0.0 P34 Ugva0 0.0 0.00 0.00 0.00 0.0 0.00 
Zi 0.0 1.34 8.40 0.0 0.00 0.00 0.00 0.0 0.00 
28 0.0 2.34 800 0.0 0.00 0.00 0.00 0.0 0.00 
29 0.0 1.34 8.40 0.0 0.00 0.00 0.00 0.0 0.00 
30 0.0 1.34 8.40 0.0 0.00 0.00 0.00 0.0 0.00 
oe 0.0 1.34 8.40 0.0 0.00 0.00 0.00 0.0 0.00 
SZ 0.0 1.34 8.40 0.0 0.00 0.00 0.00 Sage 0.00 
og 0.0 1.34 8.40 0.0 0.00 0.00 0.00 0.0 0.00 
34 0.0 1.34 . 8e740 0.0 0.00 0.00 0.00 0.0 0.00 
35 0.0 134 Seee O20 0.00 0.00 0.00 0.0 0.00 
36 0.0 1.34  VEH4a0 0.0 0.00 0.00 0.00 0.0 0500 
57 0.0 1,347 ee 80 0.0 0.00 0.00 0.00 0.0 0.00 
38 0.0 1.34 8.40 0.0 0.00 0.00 0.00 0.0 0.00 
39 0.0 1234) 98540 0.0 0.00 0.00 0.00 0.0 0.00 
40 0.0 1.34 8.40 0.0 0.00 0.00 0.00 0.0 0.00 
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ENPUT DATA FILE E2: 


SIMULATION OF FIELD EXPT. 2a: 14C-LINDANE APPLIED TO Ap HORIZON. 
WR:WMCAP SET AT 2 BARS. LAYERS 
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350 
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OUTPUT DATA FILE E2: 


PFSS2: TWO-POOL MODEL OF SOLUTE TRANSPORT THROUGH 
PARTIALLY FROZEN STRUCTURED SOIL. 


SIMULATION OF FIELD EXPT. 2a: 14C-LINDANE APPLIED TO Ap HORIZON. 
WR: WMGAP SET PAT M26 BARS © CAYERS -="40 *X"25 74MM Se SITEV TO. 

40-48: 430% S40 ee raS 82504 
(NOTE: DATA ECHO, INITIAL CONDITIONS, INFILTRATION RECORDS, 
INTERMEDIATE SOLUTE DISTRIBUTION PROFILES, ARE OMITTED) 


SOLUTE PROFILE, DAY 48 
LAYR GRAVHOH MOBHOH RETHOH GRAVSOL MOBSOL RETSOL DISSOL ADSOL TTLSOL 


1 0.0 0.0 ToS 0.0 0.0 1.49 CEN olay Aileen tsVagen 
2 0.0 0.0 Fecal’ 0.0 CFO 0.14 0.14 8.62 8.7/2 
se) 0.0 O20 i axe! 0.0 0.0 O02 O02 L358 1.6! 
4 ORM Laie oA Finders’ 0.00 OP CUMEFUTO0 0.00 0.09 0.10 
) F376 1.92 TOSS 0.00 0.00 0.00 Ooo) 02007550 .00 
6 176 192 7.38 0.00 0.00 0.00 O500>, ~0300 0.00 
if 1.00 ba ah 9.42 0.00 0.00 0.00 0700 0300 0.00 
8 One de2d 9.42 0.00 0.00 0.00 Os00% 0200 0.00 
9 0.72 Tes 9.42 0.00 0.00 0.00 0.00 0.00 0.00 
10 0.92 ie. 9.42 0.00 0.00 0.00 0200 -2%0.200 0.00 
11 P02 121 9.42 0.00 0.00 0.00 OOO 0200 0.00 
£2 1206 Lapel 9.42 0.00 0.00 0.00 G00; #70200 0.00 
13 0.91 Pate dt 9.42 0.00 0.00 0.00 C7007 ae 0200 0.00 
14 O72 eee 9.42 0.00 0.00 0.00 Cr O0 mr .00 0.00 
£5 OU72 ieee 9.42 0.00 0.00 0.00 GFOOm FeO 200). 0:00 
16 O° 72 det 9.42 0.00 0.00 0.00 GrO0. #*0200—. | 0700 
By O27 2 Bae 9.42 0.00 0.00 0.00 C700" #0700. 5.0500 
18 OFF 2 Aes og 9.42 0.00 0.00 0.00 Gz00> 0F00= 0.00 
19 Os72 AeA: 9.42 0.00 0.00 Or 00 B= O700)) 02002 7-0:,00 
20 Coie P20 9.42 0.00 0.00 0.00 0.00 ~§~0200 0.00 
21 Ox] 2 ieee 9.42 0.00 0.00 0.00 0700 0700 0.00 
Vues Oni2 I 9.42 0.00 0.00 0.00 0.00 0.00 0.00 
LS O72 ta 9.42 0.00 0.00 0.00 0.00 0.00 0.00 
24 Oz Leo 9.42 0.00 0.00 0.00 0.00 0.00 0.00 
Jie. 0,52 1.34 8.40 0.00 0.00 0.00 0.00 0.00 0.00 
26 0752 1.34 8.40 0.00 0.00 0.00 0.00 0.00 0.00 
27 O22 1 Poel ASC 0.00 0.00 0.00 0.00 0.00 0.00 
28 OFSZ 1.34 8.40 0.00 0.00 0.00 0.00 0.00 0.00 
29 One 1.34 8.40 0.00 0.00 0.00 0.00 0.00 0.00 
30 0.0 1L2OS ro sO ees 0 0.0 0.0 0.0 0.0 0.0 
ei! 0.0 Peo4 = 'So240 0.0 0.0 0.0 0.0 0.0 0.0 
32 0.0 1.34 8.40 0.0 0.0 0.0 0.0 0.0 0.0 
33 0.0 1.24 8.40 og gs, 0.0 0.0 0.0 0.0 
34 0.0 IS GRO eae 0.0 O.0 0.0 0.0 O70 0.0 
She. 0.0 L710" eee 340 0.0 0.0 0.0 0.0 0:50 0.0 
36 0.0 Ohl is. 8.40 0.0 O20 0.0 O20 0.0 0.0 
= 0.0 0.00 8.40 0.0 0.0 0.0 0.0 0.0 0.0 
38 0.0 0.00 8.40 0.0 0.0 Ose 0.0 0.0 0.0 
39 0.0 0.00 8.40 0.0 0.0 0.0 0.0 0.0 0.0 
40 0.0 Trot 8.40 0.0 0.0 0.0 O70 0.0 0.0 
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FORTRAN Program "PFSS2" 
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Ci PFSS2u=sFORTRANSPROGRAM IMPLEMENTING LAYERED 2-POOL 
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MODEL OF SOLUTE TRANSPORT THROUGH PARTIALLY 
FROZEN STRUCTURED SOIL. 


VARIABLES AND PARAMETERS: 


M, DELTAZ: NUMBER OF LAYERS, AND THICKNESS OF 
OF LAYERS USED IN SYSTEM MODELLED. 
TFIN: TIME PERIOD (DAYS) OF SIMULATION. 
Tl, T2, T3: TIMES WHEN INTERMEDIATE RESULTS CAN 
BE PRINTED. 
WR,WM,WGRAV,WT: RETAINED, MOBILE, GRAVITY AND 
TOTALS WATERS WRESPaye IN! LAYER ad. 
WMCAP, WGRCAP: MOBILE WATER CAPACITY AND GRAVITY 
WATER CAPACITY, RESP., IN LAYER, I. 
XSWMCP, XSWGCP: EXCESS CAPACITY, MOBILE WATER AND 
GRAVITY WATER, RESP., IN LAYER, I. 
MS: TOTAL MASS (ST(I)+Q(I) IN LAYER I. 
SR,SM,SGRAV,ST: SOLUTE HELD IN RETAINED, MOBILE, AND 
GRAVITY WATER POOLS, AND TOTAL DIS- 
SOLVED SOLUTE, RESP.., IN LAYER, 1. 
DB, KADS: BULK DENSITY AND ADSORPTION COEFFICIENT, 
RESPECTIVELY, USED FOR PARTITION OF 
ADSORBED SOLUTE BETWEEN ADSORBED AND 
SOLUTION PHASES. 
B: EFFECTIVE PARTITION (ADSORPTION) COEFFICIENT 
(DB*KADS/WT) . 
WA, SA: WATER AND SOLUTE TO LAYER, I, FROM ABOVE. 
WL, SL: WATER AND SOLUTE DISPLACED FROM LAYER, I, 
TO LAYER I+1, BY WA AND SA. WA AND SA ARE 
RESET EQUAL TO WL AND SL AS I IS 
INCREMENTED. 
WDRAIN, SDRAIN, DCONCN: SINK TERMS, AMT OF 
DRAINAGE WATER, AMT OF DRAINAGE SOLUTE, 
AND CONCENTRATION OF S IN WDRAIN, RESP. 
FLNR,FLNRT: FROZEN LAYER NUMBER, AND PREVIOUS 
(J-1) FROZEN LAYER NUMBER, RESP. 
FENRD ALSO,ACTS AS A ‘SWITCH’ TO 
DIRECT THE PROGRAM TO THE APPROPRIATE 
ROUTINE. 
WE, SE: WATER VOL AND SOLUTE MASS, RESP. LOST TO 
LAYER ABOVE BECAUSE OF EVAPN/CAPILLARITY 
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DOUBLE PRECISION SA, SL, DELTAS, B 

DOUBLE PRECISION WT(99), WGRAV(99), WM(99), XSWMCP (99) , 

1XSWGCP (99), SR(99), SM(99), SGRAV(99), ST(99), 

2WE (99), SE(99), MS(99), Q(99) 

REAL NETEVN, WPOND, W, E, S, WA, WL, WDRAIN, 

1DCONC, DELTAZ, DB(99), KADS(99), DEPTH(99) , 

2WGRCAP (99) , WMCAP (99), WR(99) 

DIMENSION TITLE1(80), TITLE2(80) 

WRITE (6, 608) 
608 FORMAT('PFSS2: TWO-POOL MODEL OF SOLUTE TRANSPORT THROUGH’ / 

*'PARTIALLY FROZEN STRUCTURED SOIL. '/) 

READ (S500) "TITLE: TITLEZ 

WRITE(6,500) TITLE], TITLE2 
500 FORMAT (80A1/80A1) 

READG2508) M, GERINS PY, T2003), ((DELTAZ 

WREDE (6501) “M. TFIN, T1, T2, T3, DELTAZ 
501 FORMAT (5155 F5.2) 
502 FORMAT (7F10.4) 
503 FORMAT (I3) 
504 FORMAT (3F10. 2) 
601 FORMAT(1X,'TOO MANY LAYERS -- CHECK DIMENSIONS') 
602 FORMAT('DATA ECHO: INITIAL CONDITIONS'/' DB KADS', 

xs MS WT WR WMCAP WGRCAP ' ) 
606) FORMAT (130892, 2 V2) 22% ST 202K AFG 312K . F623, 2260S, 2k, 66.3, 

*2X, F653, 20S ES) 
607 FORMAT('J FLNR FLNRT W-INPUT EVAPN S-INPUT W-POND', 

* WA(NETW) WTOTAL') 

IF (M/L7T#100)GO ‘TO: 50 

WRITE (6,601) 

GO TO 9999 
50 READ(5,502) ( (DB (RK) , KADS (K) , MS (K) , WT (K) , WR (K) , WMCAP (R) , 

*WGRCAP (K) ) ,K=1,M) 

WRITE (6, 602) 

WRITE (6,502) ( (DB (K) , KADS (K) ,MS (K) , WT (RK) , WR (K) , WMCAP (K) , 

*WGRCAP (K) ) ,K=1,M) 


Ct rr + 
C| INITIALIZE SOIL VARIABLES (WGRAV, XSWGCP, XSWMCP, WM, ST, 
G Q, SR, SM, SGRAV). 
C| INITIALIZE PROGRAM VARIABLES (J, FLNRT). 
a SS SS SS ee Se + 
WRITE (6, 2001) 
DO 51 K=1,M 


B = DELTAZ*DB (K) *KADS (K) /WT (K) 
ST(K) = MS(K)/(1+B) 

Q(K) = MS(K) *B/(1+B) 

SR(K) = WR(K) *ST(K) /WT(R) 

WM(K) = WT(K) —- WR(R) 

SM (K) WM (K) *ST (K) /WT (R) 
WGRAV(K) = 0.0 

SGRAV(K) = 0.0 

XSWGCP (K) = WGRCAP(K) - WGRAV(K) 
XSWMCP (K) WMCAP (K) - WM(R) 
DEPTH(K) = DELTAZ*(-FLOAT(K) + 0.50) 
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(MT. (on ainaeen 7 
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W, (A) AW, OOTW, OO 2M, 2) SGAA, O40) 7 (802 Spare 
iM, (2, (0) Sau 9 
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2 MWEX ,TIOWAK “ VANOW) “GseaTwaAY “oR aie 


\o “> . r 
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.(TAWIT .0) Galant dav MaRDOR? ASTaArTING 
f 008 aT aw 
“tom Fe Oa 

(AT TW Cf) aA OD) tO TATSa0 - @ 

(get) \ (A ou = ODTE 

e-|)\ae (alee = Gop 

"tw OD T2* (a) aw = Cee 

(4)5¥ - Ca) T= CA 

LO TW. CA) TE* GO Mw = OOM 

0.0 = (A) VAROW 

0.0-= (4) VARS” 

(4) VAROW~ CA)SADSOW @ On 30wae 
GOMn - (OMW (a) soewexy 

(02.0 + (OTAOJI-)*xaT isd = (ORraae 


a 


WRITE (6, 2000) (K, DEPTH (K) , XSWGCP (K) , XSWMCP (K) , WM (K) ,B, 
SOR ST CRY SRC ESMC KD 
2000 FORMAT (1X,12,F6.0,8(2X, F6.2)) 


2001 FORMAT(' I DEPTH XSWGCP XSWMCP WM B One 
a oT SR SM') 
51 CONTINUE 
eae) 
WPOND = 0.0 
FENRI= (0 


WIOTAL = 0.0 
WRITE (6,607) 
iBall a 
DO 11 K=1,M 
XSWMCP (K) = WMCAP(K) - WM(RK) 
XSWGCP (K) = WGRCAP(K) - WGRAV(K) 
11 CONTINUE 
T sol 
IF (FLNRT.GT.0)GO TO 310 
C CHECK PREVIOUS POSITION OF FROZEN LAYER AND GO TO 
C REDISTRIBUTION ROUTINE IF NECESSARY. 
READ (5,503) FLNR 
IF (FLNR.LE.M)GO TO 2 
WRITE(6,609) M, FLNR 


609 FORMAT (’SYSTEM CONTAINS ONLY ',13,'LAYERS;', 
e "SYSTEM NOT DEFINED FOR FLNR = ',13) 
GO TO 9999 


C RESET FLNRT TO FLNR AND READ IN PRECIP., EVAPN, PRECIP-BORNE 
C SOLUTE CONCENTRATION. 
2 FLNRT = FLNR 
TL=hi 
READ(5,504) W, E, S 
WIOTAL = WTOTAL + W 
WA = W - E + WPOND 
S = (S*W + SPOND*WPOND) /WA 
WRITE(6,606) J, FLNR, FLNRT, W, E, S, WPOND, WA, WTOTAL 
WPOND = 0.0 
SPOND = 0.0 
IF (WA.GE.0.0)GO TO 99 
NETEVN = —-WA 
GO TO 400 
99 IF(WGRAV(I).GT.0.0)GO TO 299 
IF(I.EQ.FLNR)GO TO 300 


(CaP SSS SSS SSIES SSO SSS SSS ar 
C| SLOW LEACHING ROUTINE. 
C+------------- - - - - - + 


100 IF(WA.GT.XSWMCP(1I))GO TO 101 
WM(I) = WM(I) + WA 
WT (I) = WR(I) + WM(I) + WGRAV(I) 
DELTAS = S*WA 
MS@) = MSU) + DELTAS 
B = DELTAZ*DB (I) *KADS (1) /WT (I) 
ST(I) = MS(I)/(1+B) 
Q(I) = MS(I)*B/(1+B) 
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C}| EQUILIBRATION OF SOLUTE AMONG RET,MOB,&GRAV POOLS: 
€| SOL GONGNGY4IS* THE SAME IN EACH WITHIN EACH LAYER. 


SR(I) = ST(1I) *WR (I) /WT (I) 
SM(1) = ST(1) *WM(I) /WT (I) 
SGRAV(I) = ST(I) *WGRAV (1) /WT (I) 


(CaF SS SES SSS SSS SSS BRS SSS SSS SS SS SS SS SS SSS + 
C RESET WA, SA, S: THESE VARIABLES CARRY WATER AND 
C SOLUTE INTO THE NEXT LAYER. 
C RESET EXCESS MOBILE WATER CAPACITY AND CONTINUE. 
Cth rrr + 
WA = 0.0 ‘ 
SAC ¢0r.0 
Spleen Cbd @, 
XSWMCP (I) = WMCAP(I) - WM(I) 
GO TO 999 


101 IF(WA.GT.WMCAP(I))GO TO 200 
WL = WA - XSWMCP (1) 
SL = ST(I) *WL/WT (1) 
SA = S*WA 
DEUBAS = SAN=TST 
MS(I) = MS(I) + DELTAS 
WM(I) = WM(I) + WA - WL 
WT(I) = WR(I) + WM(I) + WGRAV(I) 
C PARTITION OF SOLUTE, SOLUTION/ADSORBED PHASES. 
B = DELTAZ*DB (I) *KADS (1) /WT (1) 
ST(I) = MS(I)/(1+B) 
OG) | = MSGR) Bi C<B) 
C EQUILIBRATION OF DISSOLVED SOLUTE AMONG POOLS. 
SR(I) = ST(1) *WR(CI) /WT (1) 
SM(I) = ST(1I) *WM(I) /WT (I) 
SGRAV(I) = ST(1) *WGRAV (1) /WT (1) 
WA = WA - XSWMCP(1) 
SA = SL 
S = SA/WA 
XSWMCP (I) = WMCAP(I) - WM(1) 
IF(I.EQ.M)GO TO 900 


i} 


Neer i 
GO TO 99 
€ c SesGoceconece Ke teak se sevesesesesdes sesescs cise lsh) 
QRS RSE SPS SSS SS SSS SS SS SS SSS SS SS SSS SSS SS SSS + 
C| FAST LEACHING ROUTINE(200). 
(a a a a ee ar a a a ea nae oa + 
200 WL = WM(I) 
SL = SM(I) 


SA = S*WMCAP (1) 
WM (1) = WMCAP (1) 
DELTAS = SA - SL 
MS(I) = MS(I) + DELTAS 
WT(I) = WR(I) + WM(I) + WGRAV(1I) 
B = DELTAZ*DB (I) *KADS (I) /WT (1) 
C PARTITION 
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ST(I) = MS(I)/(1+B) 
Q(I) = MS(1I)*B/(1+B) 
C EQUILIBRATION 
SR(I) = ST(I) *WR (1) /WT (1) 
SM(I) = ST(1I) *WM(I) /WT (I) 
SGRAV(1) = ST(I) *WGRAV (I) /WT (1) 
C RESET WA,SA,S,XSWMCP (I) 
WA = WA — XSWMCP(1I) 
SA = SL + S*(WA-WL) 
S = SA/WA 
XSWMCP (1) = WMCAP(I) -WM(I) 
IF(I.EQ.M)GO TO 900 
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299 IF(WGRCAP (I) .GT.WGRAV(I))GO TO 301 
SOCEDLIT Et PS Td 
301 IF(I.EQ.0)GO TO 309 
IF (WA.GT.XSWGCP(I))GO TO 302 
WGRAV(I) = WGRAV(I) + WA 
DELTAS = S*WA 
MS(I) = MS(I) + DELTAS 
WT(I) = WR(I) + WM(1) + WGRAV(I) 
B = DELTAZ*DB (I) *KADS (1) /WT (1) 
C PARTITION OF SOLUTE MASS -- ADSORBED-DISSOLVED. 
ST(I) = MS(I)/(1+B) 
Q(I) = MS(I) *B/(1+B) 
C EQUILIBRATION OF DISSOLVED SOLUTE. 
SR(I) = ST(1I) *WR(I) /WT (CI) 
SM(I) = ST(1) *WM(1I) /WT (1) 
SGRAV(I) = ST(I) *WGRAV (1) /WT (I) 
C RESET FEED VARIABLES 


WA = 0.0 
SAt=2050 
5 = 000 


XSWGCP (I) = WGRCAP(I) - WGRAV(1I1) 
FLNRT = FLNR 
GO TO 999 
302 WGRAV(I) = WGRCAP (1) 
DELTAS = S*XSWGCP (1) 
MS(I) = MS(I) + DELTAS 
WT(I) = WR(I) + WM(I) + WGRAV(I) 
B = DELTAZ*DB (I) *KADS (1) /WT (1) 
C PARTITION OF SOLUTE MASS -- ADSORBED-DISSOLVED. 
ST(1) =SMS(T) /(15B) 
Q(I) = MS(I) *B/(1+B) 
C EQUILIBRATION OF DISSOLVED SOLUTE. 
SR(I) = ST(I) *WR(I) /WT (1) 
SM(I) = ST(I) *WM(I) /WT (1) 
SGRAV(I) = ST(1) *WGRAV (I) /WT (1) 
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C RESET FEED VARIABLES. 
WA = WA - XSWGCP (I) 


SA = S*WA 
S = SA/WA 
C WRITE(6,606) J, FLNR, FLNRT, W, E, S, WPOND, WA, WTOTAL 
XSWGCP (I) = WGRCAP(I) - WGRAV(I) 
GO TO 300 
Ct-- + 
C| IF WGRAV BUILDS UP ABOVE I=1, WA IS PONDED(WPOND=WA) . 
C WPOND IS THEN ADDED TO WA AT NEXT TIME STEP. 
Ct—--- 3 - nn + 
309 WPOND = WA 
SPOND = § 
WA = 0.0 
Sea0-.0 
eat 
FLNRT = FLNR 
GO TO 999 
Ct------- A A a a te 


C REDISTRIBUTION OF WATER AND SOLUTE WHEN FROST LAYER 
c MOVES DOWNWARD BETWEEN SUCCESSIVE TIME STEPS. 


310 READ(5,503) FLNR 
C IF FROZEN LAYER HAS MELTED, CLEAR WGRAV FROM THE PROFILE. 
IF(FLNR.EQ.0)GO TO 314 
C IF FROZEN LAYER HAS NOT MOVED DOWNWARD THERE IS NO REDISTRIBUTN. 
IF(FLNR.LE.FLNRT)GO TO 2 
C IF THERE IS NO WGRAV TO REDISTRIBUTE RETURN TO 2. 
DO 312 K=1,FLNRT 
IF (WGRAV (K) .GT.0.0)GO TO 314 
312 CONTINUE 
GO TO 2 
C REDISTRIBUTE WGRAV, STARTING AT THE LAYER ABOVE THE PREVIOUS FROZEN LAYER 
314 IR = FLNRT - 1 
IF(IR.EQ.0)GO TO 2 
WA = WGRAV (IR) 
IF (WGRAV (IR) .EQ.0.0)GO TO 320 
S = SGRAV (IR) /WGRAV (IR) 
WGRAV(IR) = 0.0 
XSWGCP (IR) = WGRCAP(IR) - WGRAV(IR) 
WT(IR) = WRCIR) + WMCIR) 
MSCIR) = MSCIR) - SGRAVCIR) 
B = DELTAZ*DB (IR) *KADS (IR) /WT CIR) 
C PARTITION 
ST(IR) = MS(IR)/(1+B) 
QCIR) = MSCIR) *B/(1+B) 
C EQUILIBRATION 
SR(IR) = WRCIR) *ST (IR) /WT CIR) 
SM(IR) = WMCIR) *ST CIR) /WT CIR) 
SGRAV (IR) = WGRAV (CIR) *ST CIR) /WT (IR) 
= IRS 
FLNRT = -1l 
GO TO 99 
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400 


4001 


4002 


474 
475 


401 


402 


te ate sto ale ale ale ote slo ole ol. 
TIE FE FE FE FE FEE ES 


E TO THE NEXT HIGHER LAYER AND REDISTRIBUTE WGRAV DOWNWARD. 


IR = IR - 1 
TE GIRGEQ.0)) GO TO 2 
WA = WGRAV (IR) 
IF (WGRAV (IR) .EQ.0.0)GO TO 320 
S = SGRAV (IR) /WGRAV (IR) 
WGRAV(IR) = 0.0 
XSWGCP (IR) = WGRCAP(IR) - WGRAVCIR) 
WT(IR) = WRCIR) + WMCIR) 
MS(IR) = MS(IR) - SGRAV(CIR) 
B = DELTAZ*DB (IR) *KADS (IR) /WT CIR) 
ST(IR) = MSCIR) /(1+B) 
QCIR) = MS(CIR) *B/(1+B) 
SR(IR) = WRCIR) *ST CIR) /WT CIR) 
SM(IR) = WM(IR) *ST (CIR) /WT (IR) 
SGRAV(IR) = WGRAV (CIR) *ST (CIR) /WT CIR) 
1 *= BERL /1 
GO TO 99 


0 


VAPORATION ROUTINE: WATER IS REMOVED FROM 

TOP LAYER (I=1) FIRST. WATER AND SOLUTE 

MOVE UP THE PROFILE TO REPLENISH WR(1) ONLY. 
SOLUTE VOLATILIZED FROM LAYER 1 ONLY. 

UPWARD MOVEMENT CONTINUES UNTIL NET EVAPORATN 
LOSSES ARE SATISFIED BY WGRAV + WM. 

UPWARD MOVEMENT OF WATER IS ACCOMPANIED BY 
UPWARD MOVEMENT OF SOLUTE. 

WHERE NETEVN IS NOT SATISFIED BY WT(1), THE 
PROCESSES ABOVE ARE REPEATED UNTIL NETEVN = O. 


IF(NETEVN .LE .WT(1)) GO TO 4001 
NETEVN = NETEVN - WT(1) 
WE(1) = WT(1) 
GO ‘TO, "4002 
WE(1) = NETEVN 
NETEVN = 0.0 
SE(1) = WE(1)*ST(1)/WT(1) 
WRITE (6, 474) 
FORMAT('NET E WT(1) WG(1) WM(1) WR(1) WE(1)~ SE(1)') 
FORMAT (7 (F5.2, 2X) ) 
WRITE (6, 475) NETEVN,WT (1) ,WGRAV (1) ,WM(1) ,WR(1) ,WE(1) , SE(I) 
WGMSUM = WGRAV(1) + WM(1) 
I=1l 
IF(WE(1) .LE. WGMSUM)GO TO 402 
De + S 
WGMSUM = WGMSUM + WGRAV(I) + WM(1I) 
CO ENG 407 
ELLAST = I 
IF (ELLAST.EQ.1)GO TO 405 
K = ELLAST - 1 
DO 403 L=1,K 
LNEXT = L + 1 
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WE (LNEXT) 
SE (LNEXT) 
403 CONTINUE 
K = ELLAST - 1 
DO 404 L=1,K 
LNEXT. =@ 4 
WT(L) = WT(L) - WE(L) + WE(LNEXT) 
DELTAS = -SE(L) + SE(LNEXT) 
MS(L) = MS(L) + DELTAS 
WM(L) = 0.0 
SM(L) = 0.0 
WGRAV(L) = 0.0 
SGRAV(L) = 0.0 
WR(L) = WT(L) 
B = DELTAZ*DB (L) *KADS (L) /WT (L) 
PARTITION OF SOLUTE MASS. 
ST(L) = MS(L) /(1+B) 
OCE) = MSL) «By (1t8) 
SR(L) = ST() 
404 CONTINUE 
405 WT(ELLAST) = WT(ELLAST) - WE(ELLAST) 
DELTAS = -SE(ELLAST) 
MS (ELLAST) = MS(ELLAST) + DELTAS 
B = DELTAZ*DB (ELLAST) *KADS (ELLAST) /WT (ELLAST) 
PARTITION OF SOLUTE MASS. 
ST(ELLAST) = MS(ELLAST) / (1+B) 
Q(ELLAST) = MS(ELLAST) *B/ (1+B) 
WGRAV(ELLAST) = WT(ELLAST) - WM(ELLAST) —- WR(ELLAST) 
IF (WGRAV (ELLAST) .LT.0.0)GO TO 406 
WM (ELLAST) = WMCAP (ELLAST) 
GO TO 499 
406 WGRAV(ELLAST) = 0.0 
WM(ELLAST) = WT(ELLAST) - WR(ELLAST) 
499 SGRAV(ELLAST) = WGRAV(ELLAST) *ST (ELLAST) /WT (ELLAST) 
SM(ELLAST) = WM(ELLAST) *ST(ELLAST) /WT (ELLAST) 
SR(ELLAST) = WR(ELLAST) *ST(ELLAST) /WT (ELLAST) 
IF(NETEVN .GT. 0.0)GO TO 400 
GO TO 999 
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WE(L) - WM(L) - WGRAV(L) 
WE (LNEXT) *ST (LNEXT) /WT (LNEXT) 


DRAINAGE ROUTINE: IF WATER IS LOST FROM THE BOTTOM LAYER, 
IT IS REMOVED FROM THE SYSTEM AS DRAINAGE. 


900 WDRAIN = WA 
SDRAIN 
DCONCN 
WRITE (6,603) J 
WRITE (6,604) 

WRITE(6,605) WDRAIN, SDRAIN, DCONCN 

603 FORMAT(2X, ‘DRAINAGE WATERS AFTER DAY',I2) 

604 FORMAT (' VOLUME AMT SOLUTE CONC. SOLUTE,’) 

605 FORMAT (3F10.3, 2X) 

999 IF(IR.GT.0)GO TO 320 


now 
nw 
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699 
698 


697 


SESS, 


IF (FLNRT.LT.0)GO TO 2 

FLNRT = FLNR 

IF (J.NE.T1.AND.J.NE.T2.AND.J.NE.T3.AND.J.NE.TFIN)GO TO 1 

WRITE (6,699) J 

WRITE (6, 698) 

WRITE (6,697) ((K, WGRAV (K) , WM (K) , WR (K) , SGRAV (K) , SM(K) , SR(K) , ST(K 


=), OKT MSCK) ey Kad, M) 


FORMAT (1X/,'SOLUTE PROFILE, DAY ',I2) 
FORMAT ('LAYR GRAVHOH MOBHOH RETHOH GRAVSOL MOBSOL' 


oR EO SOLmeDISSOL. sADSOL! TILSOL ) 


FORMAT (TXou2s 2anhor2, 24, FO. 2,26, bOsey eae Ool sy 2m FOe2 2k, EOw2 


RD KGEG UDO Xe TOD ek Fore) 


WRITE (6,607) 
IF(J.LT.TFIN)GO TO 1 
STOP 

END 
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(sr, “oh eg ' rat (xt) 
‘102d0M J02VAND. HONTEA ZOWEOH DAYAL’ 


('SO847T JoRCA — i 
£47, XS, 0,87, 48, 5.99, 48, 5.99, MSS. BF 28,2 Bao, 1k ! 
(£.04, 48 Sabah 
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